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1.0  SUMMARY 

This  report  contains  the  development  of  an  unstructured  gnd— based  finite— volume 
integration  scheme  for  solving  the  time-domain  Maxwell’s  equations  to  study  a  myriad  of 
problems  in  electromagnetics.  The  principal  application  of  this  technology  is  the  prediction 
of  radar  cross  section  (RCS)  of  low  observable  platforms.  This  work  was  performed  under 
the  AFOSR  contract  F49620-93-C-0066. 
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1.0  INTRODUCTION 


The  design  and  development  of  modern  aerospace  configurations  increasingly  require 
better  understanding,  management  and  exploitation  of  ever  more  complex  physical  phe¬ 
nomena  in  different  disciplines  such  eis  fluid  dynamics,  structural  mechanics,  propulsion, 
controls,  and  electromagnetics  that  all  play  an  interdisciplinary  role.  The  traditional  ap¬ 
proach  of  serial  design  by  individual  disciplines  which  cannot  account  for  coupling  effects 
will  not  only  result  in  poor  design  but  also  is  not  cost  effective.  Extensive  simulations 
of  the  veuious  complex  phenomena  are  required  to  imderstand  the  interdisciplinary  role 
in  design.  Along  with  advances  in  experimental  simulations,  with  the  emergence  of  the 
supercomputers,  both  conventional  (vector  Cray-like  architectures)  and  massively  parallel, 
the  computing  technology  is  beginning  to  play  an  ever  increasing  role  in  design  simulations. 
Both  government  and  industry  view  ‘Computational  Sciences’,  a  discipline  that  exploits 
advances  in  numerical  algorithms  development  and  the  increasing  power  of  supercomput¬ 
ers,  super-minicomputers  and  graphics  workstations,  as  a  critical,  potentially  efficient  and 
cost-effective  technology  for  advanced  design. 

The  development  of  a  computational  environment  that  encompasses  different  dis¬ 
ciplines  for  multidisciplinary  studies  as  described  in  Fig.  1  requires  that  the  following 
computational  capabilities  be  properly  matured  within  each  discipline. 

1)  Computational  Fluid  Dynamics  (CFD) 

—  Turbulence  modeling 
—  Transition 
—  Finite-rate  chemistry 

—  Algorithms  for  incompressible  to  hypersonic  Mach  number  range 
—  Unsteady  and  separated  flows 

2)  Computational  Electromagnetics  (CEM) 

—  Proper  implementation  of  Miixwell’s  equations  with  general  material  properties 
—  Formulation  and  implementation  of  different  boundary  conditions 

3)  Computational  Structural  Mechanics  (CSM) 

—  Accurate  finite-element  models  for  static  and  dynamic  flexible  effects 
—  Failure  and  fatigue  analysis  accounting  for  material  grain  structure 

The  development  of  computational  capabilities  in  all  these  disciplines  critically  de¬ 
pends  on  the  following  technologies. 

1)  Geometry  and  Grid  Setup 

—  Three-dimensional  geometry  modeling 
—  Structured  and  unstructured  grid  cell  representation 
—  Pre  and  post  processing  for  visualization 
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2)  Computer  Architectures 

—  Vector /parallel  coarse  grain  machines  such  as  the  Cray  C-90  with  16  CPUs 
—  Massively  parallel  machines,  both  SIMD  (Single  Instruction  Multiple  Data)  and 
MIMD  (Multiple  Instruction  Multiple  Data) 

—  Advanced  graphics  workstations 

One  of  the  computational  disciplines  that  has  been  a  supercomputing  pace  setter  for 
the  last  three  decades  is  CFD.  This  technology,  which  started  with  the  development  of 
the  transonic  small- disturbance  theory  in  the  late  60’s,  has  matured  both  in  algorithm 
and  code  development  to  the  point  of  today  being  able  to  solve  the  time  averaged  Navier- 
Stokes  equations  for  predicting  the  flowfield  over  a  complete  fighter.  Today,  CFD  is  playing 
a  critical  role  in  the  development  of  next  generation  fighters  and  the  National  Aerospace 
Plane,  though  only  in  single  discipline  mode.  The  development  of  a  multidisciplinary 
computational  environment  can  significantly  benefit  by  applying  many  of  the  attributes  of 
CFD  to  other  disciplines  such  as  CEM. 

1.1  Attributes  of  CFD 

1)  The  fluid  dynamic  equations  are  usually  cast  in  conservation  form  either  as  differen¬ 
tial  or  local  integral  equations  conserving  mass,  momentum,  and  energy  fluxes,  thus 
allowing  for  numerical  capture  of  flow  discontinuities  such  as  shocks  and  slip  surfaces. 
Equations  representing  the  physics  in  other  disciplines,  for  example  Maxwell’s  equa¬ 
tions  in  electromagnetics,  can  be  cast  in  similar  conservation  forms  for  conservation 
of  appropiate  fluxes. 

2)  Recent  developments  of  hyperbolic  algorithms  for  solving  the  time- domain  Euler  equa¬ 
tions  are  based  on  the  characteristic  theory  of  signal  propagation  and  are  referred  to 
as  the  ‘upwind’  schemes.  For  hyperbolic  equations,  the  upwind  based  schemes  can 
be  constructed  to  provide  the  right  amount  of  numerical  dissipation  to  achieve  both 
stability  and  accuracy.  Current  state  of  the  art  in  numerical  algorithms  is  based  on 
Essentially  Non  Oscillatory  (ENO)  interpolation  schemes  that  can  provide  arbitrarily 
high  order  of  accuracy  for  arbitrary  cell  shapes  such  as  hexahedral,  triangular  prism, 
and  tetrahedral  elements. 

3)  For  treatment  of  complex  aerospace  configurations,  CFD  methods  usually  employ 
either  a  structured  grid  based  body-fitted  coordinate  system  or  an  unstructured  finite- 
element  grid  setup  for  ease  in  implementing  the  various  boundary  conditions.  Similar 
numerical  geometry  and  grid  setup  procedures  are  equally  applicable  to  modeling 
complex  problems  in  other  disciplines. 

4)  Pre  and  post  processing  capabilities  running  on  advanced  graphics  work  stations  are 
effectively  employed  to  visualize  and  animate  the  geometry/grid  and  solution. 

The  goal  of  ‘Computational  Sciences’  is  to  effectively  employ  the  many  advances  in 
CFD  coupled  with  emerging  supercomputer  architectures  with  expected  teraflops  (trillion 
floating  point  operations  per  second)  performance  to  mature  the  computational  technology 
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in  different  disciplines  and  be  able  to  perform  multidisciplinary  studies  critical  to  advanced 
design.  " 

1.2  Objectives 

Toward  establishing  a  computational  environment  for  performing  multidisciplinary 
studies,  the  initial  goal  is  to  advance  the  state-of-the-art  in  CEM  with  the  following  specific 
objectives. 

1)  Apply  algorithmic  advances  in  CFD  to  solve  Maxwell’s  equations  in  general  form  to 
study  scattering  (radar  cross  section),  radiation  (antenna),  and  a  variety  of  eletro- 
magnetic  environmental  (electromagnetic  compatibility,  shielding,  and  interference) 
problems  of  interest  to  both  the  defense  and  commercial  community. 

2)  Mature  the  CEM  technology  to  the  point  of  being  able  to  perform  coupled  CFD/CEM 
optimization  design  studies. 

3)  Establish  the  viability  of  MIMD  massively  parallel  architectures  for  tackling  large 
scale  problems  not  amenable  to  present  day  supercomputers. 

1.3  Computational  Electromagnetics  (CEM) 


The  ability  to  predict  radar  return  from  complex  structures  with  layered  material 
media  over  a  wide  frequency  range  (100  MHz  to  20  GHz)  is  a  critical  technology  need 
for  the  development  of  stealth  aerospace  configurations.  Traditionally,  radar  cross  section 
(RCS)  calculations  have  employed  one  of  two  methods:  high  frequency  eisymptotics,  which 
treats  scattering  and  diffraction  as  local  phenomena;  or  solution  of  an  integral  equation 
(in  the  frequency  domain)  for  radiating  sources  on  (or  inside)  the  scattering  body,  which 
couples  all  parts  of  the  body  through  a  multiple  scattering  process.  A  third  approach 
is  the  direct  integration  of  the  differential  or  integral  form  of  Maxwell’s  equations  in  the 
time-domain. 

The  time-domain  Maxwell’s  equations  represent  a  more  general  form  than  the 
frequency-domain  vector  Helmholtz  equations,  which  are  usually  employed  in  solving  scat¬ 
tering  problems.  A  time-domain  approach  can,  for  instance,  handle  continuous  wave 
(single  frequency)  as  well  as  a  single  pulse  (broadband  frequency)  transient  response. 
Frequency-domain-based  methods  usually  provide  the  RCS  response  for  all  angles  of  inci¬ 
dence  at  a  single  frequency,  while  time-domain-based  methods  provide  solutions  for  many 
frequencies  from  a  single  transient  calculation.  Also,  in  a  time-domain  approach,  one 
can  consider  time-varying  material  properties  for  treatment  of  active  surfaces.  By  using 
Fourier  transforms,  the  time-domain  transient  solutions  can  be  processed  to  provide  the 
frequency-domain  response.  Frequency-dependent  (dispersive  )  and  anisotropic  material 
properties  can  also  be  included  within  the  time-domain  formulation. 

CEM  is  a  critical  technology  for  the  United  States  to  maintain  its  leadership  in  the 
advancement  of  future  aerospace  development  through  supercomputing.  As  we  transition 
from  the  present  Gigaflops  to  the  next  generation  Teraflops  computing,  CEM  will  become 
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integral  to  aerospace  design  not  only  as  a  stand  alone  technology  but  also  as  part  of  the 
multidisciplinary  coupling  that  leads  to  well  optimized  designs. 


1.4  CEM  Issues 

Proper  development  of  a  CEM  capability  appropriate  for  ^1  aspects  of  aerospace 
design  must  consider  various  issues  associated  with  electromagnetics.  Some  of  them  are. 

1.  Physics  of  Maxwell’s  equations 

•  differential  and  integral  forms  of  Maxwell’s  equations  suitable  for  numerically  satisfy¬ 
ing  tangential  field  flux  conservation 

•  flexibility  in  implementing  total  fleld,  scattered  field,  and  diffracted  field  forms  de¬ 
pending  on  the  nature  of  the  problem  being  solved 

•  incorporation  of  various  source  terms 

2.  Material  properties 

•  perfectly  conducting  surfaces 

•  lossy/lossless  e  and 

•  resistive  sheet  (conductivity  (t  and  thickness  d) 

•  impedance  layer 

•  anisotropic  media 

•  chiral  media 

•  dispersive  media  -  e(u))  and 

•  nonlinear  materials 

3.  Boundary  conditions 

•  perfectly  conducting,  nxS  =  0  (Electric  wall)  and  n  x  W  =  0  (magnetic  wall) 

—  accurate  evaluation  of  n  x  7Y  on  the  PC  surface  is  crucial 

•  material  interface,  |n  x  S\  and  |n  x  'H\  are  zero 

—  algorithms  must  account  for  any  variation  in  e  and  ^  at  interface 

•  resistive  sheet,  n  x  S  and  n  x  Ti  jump  across  RS 

•  impedance  layer,  n  x  n  x  £  =  —rj  nxH 

—  implementation  in  time— domain  involves  convolution  integrals 

•  nonreflecting  farfield 

—  characteristics  based  hierarchy  of  absorbing  conditions  for  total  field,  scattered 
field,  and  diffracted  field  formulations 

•  periodic 

•  zero  flux  (collapsing  cell  faces) 

4.  Algorithmic  issues 


5 


•  unstructured  grid-based  finite-volume  algorithms  that  include  structured  grids  as  a 
speaal  case 

•  stability  and  accuracy  of  schemes  using  spectral  techniques 

•  construction  of  higher— order  schemes  including  at  boundauies  using  polynomial  rep¬ 
resentations  for  electric  and  magnetic  field  vziriations  inside  a  general  polyhedron  cell 

•  implicit  and  explicit  schemes  with  approprite  space  and  time  discretization 

5.  Gridding 

•  geometry  modeling  using  CAD  packages 

•  structured  or  unstructured  surface  grids  using  advancing  front  techmque 

•  volume  grids  using  hexahedreil,  prismatic  or  tetrahedral  cells 

•  optimization  routines  for  achieving  grid  quality  and  smoothness 

•  adaptive  gridding 

6.  Massively  parallel  computing  issues 

•  domain  decomposition  and  load  balancing 

•  internodal  message  passing  with  minimum  communication  delays 

•  synchronization  for  time  accurate  computation 

•  measure  of  MFLOP  rating 

•  scalability  measures 

•  pre  and  post  processing  in  parallel  environment 

•  transportability 

7.  User  Issues 

•  problem  set  up  and  boundary  conditions  -  how  flexible  and  general  is  the  code? 

•  geometry  and  grid  set  up  time,  resolution  requirements  and  computational  domain 
size 

•  intern2il  consistency  check  for  spotting  user  errors  and  dignostic  measures 

•  run  enviroment 

—  selection  of  number  of  nodes 

—  load  balancing  and  domain  decomposition 

—  automatic  termination  criteria 

—  complete  monostatic  runs 

—  post  processing  routines  -  FFT,  plotting,-  •  • 

•  Reliability  of  solution 

8.  Validation  and  Applications 
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code  validation  on  Electromagnetic  Code  Consortium  test  cases  and  canonical  solu¬ 
tions 

radar  cross  section  of  low  observable  platforms 
antenna  performance 

Microwave  monolithic  integrated  circuit  (MMIC)  modeling  and  photonic  band  gap 
periodic  structures 

electromagnetic  environmental  effects  (E®),  such  as  EMP,  EMI,  and  compatibility 
problems 

bioelectromagnetics  such  as  microwave  hyperthermia  cancer  treatment  of  humans  and 
effects  of  cellular  phones 

multidisciplinary  problems  invloving  optimization  of  performance  for  electromagnet¬ 
ics  (stealth),  aerodynamics  (flow  management),  structures  (fatigue  and  failure),  and 
controls 


2.  MAXWELL^S  EQUATIONS 

A  general  theoretical  framework  for  electromagnetic  scattering  in  the  time  domain  is 
provided  by  Maxwell’s  equations  relating  the  spatial  derivatives  of  the  electric  and  magnetic 
fields  to  their  time  derivatives  and  to  both  external  and  internal  sources.  In  conventional 
SI  notation,  the  Maxwell  curl  equations  in  the  presence  of  polarizable  materials  are  written 
as: 


VxB+^=0  , 

(1) 

(2) 

where  E  and  H  are  the  electric  and  magnetic  field  intensities,  respectively,  while  D  and 
B  are  the  electric  displacement  and  magnetic  induction,  and  J  is  the  electric  current  of 
“free”  charges  (to  be  defined  below).  These  equations  axe  supplemented  by  two  divergence 
conditions: 

V  •  £  =  0  ,  and  (3) 

V-D  =  p  ,  (4) 

where  p  is  the  volume  density  of  “free”  electric  charges.  Taking  the  divergence  of  Eq.  (1), 
one  finds  that  5(V  •  B)/dt  =  0,  so  that  Eq.  (3)  can  be  treated  as  an  initial  condition  on 
the  fields  that  will  automatically  hold  at  all  later  times. 

Similarly,  the  divergence  of  Eq.  (2)  gives: 

'^-J  +  d{V-D)ldt  =  0  ,  (5) 

while  conservation  of  the  “free”  electric  charge  requires  that 

V  ■  J  +  dp/dt  =  0  .  (6) 

Combining  Eq.  (5)  and  Eq.  (6)  gives: 

d{V  -D-  p)ldt  =  0  ,  (7) 

so  that  iiV-D  —  p&t  some  initial  instant,  Eq.  (4)  will  also  hold  at  all  future  time. 

Each  polarizable  material  is  characterized  by  a  volume  density  of  electric  dipole  m<> 
ment  P  and  magnetic  dipole  moment  M  in  terms  of  which  the  electric  displacement  D 
and  the  magnetic  field  intensity  H  can  be  written  as: 

D  =  to  E  P 

-  1  -  - 
H  =  -B-M 

po 
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(8) 

(9) 


where  eo  is  the  permittivity  of  free  space  8.854  x  10  farad/met^  in  SI  umts),  and /io 
is  the  penneabiUty  of  free  space  (  =  47r  x  10"'^  henry/meter  in  SI  units).  The  polarization 
P  and  M  in  general  can  arise  from  many  types  of  forces  on  the  material  medium,  but  in 
radar  scattering  the  most  important  of  these  is  a  linear  response  to  the  local  electric  field 
E  and  magnetic  induction  B.  The  same  is  true  of  the  current  of  “free”  charges  J  induced 
by  these  fields. 

In  this  limit;  the  most  general  form  for  the  response  R  (standing  for  either  P,  M,  or 
J)  can  be  written  as; 

R{f,t)  =  R{F,0)+  [  [r.(r,t-T)E{r,r)  +  n(f,t-T)B(,f,r)]dT  .  (10) 

Here,  the  properties  of  the  material  are  lumped  intojhe  susceptibilities  r*  and  r*.  Each 
Ta  is’ a  tensor,  reflecting  the  fact  that  the  response  R  need  not  be  parallel  to  the  applied 
field.  Note  also  that  each  ra  depends  only  on  the  time  difference,  t  -  r,  so  that  the 
integral  in  equation  (10)  is  a  simple  convolution.  The  term  P(f,0)  includes  any  static, 
field-independent  polarization  or  current,  as  one  finds  in  permanent  magnets  or  supercon¬ 
ductors,  as  well  as  any  polarizations  due  to  fields  present  before  t  =  0. 

These  additional  relations,  together  with  the  definitions  D  =  eoE  +  P  Bud  H  = 
-^B  -  constitute  a  complete  linear  system  of  equations  that  can  be  solved  for  the 

development  of  E  and  P  from  some  initial  state  that  satisfies  the  subsidiary  conditions 
V-.B  =  0andV'.D  =  p/.  Equations'(8),  (9),  and  (10)  together  already  contain  implicitly 
the  boundary  conditions  on  the  fields  that  must  be  satisfied,  for  instance,  at  the  interface 
between  material  (P  and  M  nonzero)  and  vacuum  (P  M  0). 

In  many  cases  of  interest,  the  time  scale  for  the  development  of  E  and  P  is  much 
longer  than  the  time  it  takes  for  the  material  to  build  up  its  response.  In  this  limit,  the 
relations  (10)  become  approximately  instantaneous,  and  one  can  write; 

P  =  PtE  PmB  , 

M  =  meE  +  mmB  , 

J  =  jcE  +  jmB  , 

which  greatly  simplifies  the  solution  process.  For  ordinary  materials,  a  further  simpli¬ 
fication  comes  from  the  fact  that  a  magnetic  field  does  not  produce  either  an  electric 
polarization  {pm  =  0)  or  a  local  current  {j„,  =  0),  nor  does  an  electric  field  induce  any 
magnetization  (me  =  0).  (Materials  for  which  pm  and  me  nonzero  are  called  chiral.) 
This  allows  one  to  write  directly  simple  expressions  for  £>,  P,  and  J; 

D  =  eoE+PeE  =  eE  , 

H  = —B -rrimB  =  p~^B  ,  (^2) 

po 

J  =  aE  , 
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in  terms  of  new  tensor  material  parameters  e  (the  permittivity),  /i  (the  permeability),  and 
a  (the  conductivity),  eliminating  the  need  to  calculate  P  and  M  separately.  Equation  (2) 
then  becomes: 

V  X  (m-'p)  =  ^  ,  (13) 

while  equation  (1)  is  unchanged. 

Because  equations  (1),  (2),  and  (10)  are  linear,  one  can  Fourier  transform  these  equa¬ 
tions  with  respect  to  time,  malting  use  of  the  convolution  theorem  to  obtain  a  form  similar 
to  (11),  (12),  and  (13): 

P{u})  -  pe(u;)E(u)  ;  D(u;)  =  e{u})E{u}) 

M{u)  =  fhm{(^)B{u})  ;  H{u)  =  (14) 

V  X  P  =  iuB  ;  V  X  P  =  -weE  -1-  ^E  , 

where  we  have  written  A{tjj)  =  for  every  quantity  in  (1),  (2),  and  (10). 

This  is  the  standard  form  in  which  the  material  parameters  e,  /x,  and  a  are  measured  and 
reported,  using  applied  fields  of  a  well  defined  frequency  w. 
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3.  GENERAL  CONSERVATION-LAW  FORM 


For  many  problems  in  mathematical  physics,  the  physical  process  to  be  modeled  is 
governed  by  an  appropriate  set  of  linear  or  nonlinear  partial  differential  equations.  For 
example,  many  fluid  dynamic  processes  are  governed  by  the  Navier-Stokes  equations,  and 
the  electromagnetic  scattering  from  objects  is  modeled  by  Maxwell’s  equations. 

In  general,  many  of  these  equations  naturally  lend  themselves  to  a  conservation  form 
representation  given  by: 

Qt  +  +  .Fy  +  =  ^(Source)  ,  (15) 

where  the  dependent  variable  vector  Q,  and  the  fluxes  E,  T ,  and  Q  and  the  source  S  take 
on  different  forms  depending  on  the  physical  process  being  modeled.  The  integral  form 
of  the  conservation  laws  can  easily  be  derived  from  the  differential  form  by  integrating 
Eq.  (15)  with  respect  to  x,y,z  over  any  conservation  cell  whose  volume  is  V: 


This  can  be  rewritten  in  vector  notation  as: 


|///^ec(x<iyfc+///^(v.F) 


dx  dydz  =  S 


(17) 


In  the  above: 

F  =  £j+J^k  +  Ql 


(18) 


where  j,  fc,  and  /  are  unit  vectors  along  the  i,  y,  and  z  directions,  respectively. 

Applying  the  Gauss  divergence  theorem,  we  can  convert  the  volume  integral  into  an 
integral  over  the  surface  A  that  bounds  the  volume: 


|(QV)+//^(#n)iA  =  5  .  (19) 


In  the  above  equation,  the  cell  average  of  the  dependent  variables  is  denoted  by  Q.  The 
outward  unit  normal  at  any  point  of  the  boundary  surface  of  a  cell  has  been  denoted  by 


n 


fixj  "h  Tiyk  “h 


(20) 


The  integral  form  of  the  conservation  laws  given  by  Eq.  (19)  defines  a  system  of  equations 
for  the  cell  average  values  of  the  dependent  variables.  In  order  to  construct  numerical 
methods  to  solve  the  integral  form  of  the  conservation  laws,  one  must  be  able  to  define  cell 
geometries,  approximate  the  dependent  variables,  develop  spatial  discretization  procedures, 
and  develop  time  integration  procedures  to  update  the  cell  averages,  etc.  There  are  many 
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numerical  algorithmic  issues  that  come  into  play  in  devising  a  solution  procedure  to  solve 
either  the  differential  form,  Eq.  (16),  or  the  integral  form,  Eq.  (19).  Some  of  them  axe 
1)  implicit  and  explicit  schemes,  2)  stability  and  order  of  accuracy,  3)  relaxation  and 
approximate  factorization  procedures,  4)  central  differenced  and  upwind  schemes,  5)  finite- 
volume  and  finite  difference  schemes  applied  to  a  structured  grid  (usually  for  the  differential 
form),  and  6)  finite-element-like  finite-volume  schemes  for  an  unstructured  grid  (applied 
to  the  integral  form)  setup. 

Application  of  Eq.  (15)  to  many  realistic  problems  requires  a  coordinate  transforma¬ 
tion  to  properly  represent  the  physical  domain  of  interest  and  to  aid  in  the  boundary 
condition  treatment. 

Under  the  transformation  of  coordinates  implied  by: 

T  =  t,  ^  =  f]  =  'n{t,x,y,z),  C  =  C{t,x,y,z)  , 


(21) 


(22) 


(23) 

and  6,  (y,  rjt,  J?!,  rjy,  rj,,  0,  Cx,  and  G  are  the  transformation  metrics. 

Associating  the  subscripts  j,k,l  with  the  G  C  directions,  a  numerical  approximation 
to  Eq.  (21)  (as  well  as  Eq.  (19))  may  be  expressed  in  the  semidiscrete  conservation  law 
form  given  by: 

+  (^j,k+i/2,l  -  (24) 

+  (pj.k, 1+112  -  Qj,k,l-l/2^  =  S  , 


Eq.  (15)  can  be  recast  in  the  conservation  form  given  by: 

Qr  +  ~ 5 

where:  ^ 

—  Q 

: 

7=^Q  +  '^£+’^J^+jQ  ■ 
V=jQ  +  j£  +  ^y^+jS  . 

where,  in  turn,  J  is  the  Jacobian  of  the  transformation: 

J  =  d{T,^,T],C)/d{t,x,y,z)  =  ^ 
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where  6,  andQ  are  the  numerical  fluxes  representing  the  physical  fluxes  5,  T ,  and  Q  at 
the  bounding  sides  of  the  cell  for  which  discrete  conservation  is  considered,  and  Qi.fc,/  is 
the  cell  average  of  the  dependent  variables.  The  half-integer  subscripts  denote  cell  sides 
and  the  integer  subscripts  the  cell  itself  or  its  centroid. 

The  semidiscrete  conservation  law  form  given  by  Eq.  (24)  may  be  regarded  as  rep¬ 
resenting  a  finite  volume  discretization  because  the  Jacobian  J  appearing  in  Eq.  (23)  is 
associated  ■with  the  volume  of  the  cell,  and  the  metrics  and  so  on  appearing  in  the 

fiux  terms  (Eq.  (22))  axe  nothing  but  the  components  of  the  appropriate  cell  surface  area. 

The  objective  is  to  solve  Eq.  (24)  for  the  dependent  vector  Q.  After  incorporation  of 
proper  fiux  representation,  the  discrete  form  of  Eq.  (24)  can  be  written  as: 

iJ(Q)  =  0  .  (25) 


If  Q  is  sought  in  the  neighborhood  of  a  known  solution  Q*,  then  solution  to  Eq.  (25)  can 
be  written  as:  ^ 

where  |g,  in  general,  is  a  differential  operator.  All  of  the  above-mentioned  algorithmic 
issues  apply  to  how  one  models  the  operator.  References  1-29  provide  many  of  the 
algorithmic  details. 

Unstructured  Grid  Approach 

In  the  unstructured  grid  approach,  there  is  no  defined  set  of  coordinates,  such  as  the 
J7,  C  system  in  the  differential  approach.  The  cell  shape  is  arbitrary  and  can  be  made  up 
of  any  polyhedron  shape  such  as  a  hexahedron,  prism,  or  a  tetrahedron.  The  coordinate 
directions  for  each  cell  are  defined  by  its  cell  interface  normals,  and  one  will  employ  the 
integral  conservation  form  of  equations,  Eq.  (19),  to  enforce  the  flux  conservation. 

We  first  divide  the  surface  integreil  into  component  parts  that  apply  over  each  distinct 
face  or  side  of  any  cell  under  consideration.  We  then  replace  the  integral  on  each  face  with 
a  numerical  quadrature  sis  part  of  the  numerical  approximation  procedure. 

(F  •  n)  ^  j  J  {F  ■n)dS 

faces  ^  (26a) 

=  E  E 

faces  quads. 


Here,  “quads.”  is  an  abbreviation  for  “quadrature  points.”  The  weights  of  the  quadrature 
formulae  must  include  the  effect  of  cell  surface  area  corresponding  to  the  given  face  and 
such  weights  are  denoted  by  Si  in  the  above  equation.  The  midpoint  formula  can  be 
represented  as 

f  [(F-n)dS  =  Y,  Fm-fimSm  (26&) 

^  faces 
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where  m  denotes  the  centroid  of  each  face.  Even  higher-order  quadrature  formulae  may 
be  used  if  necessary. 

The  original  initial  value  problem  (IVP)  for  the  differential  form  of  the  conservation 
laws  specifies  initial  values  of  the  dependent  variables.  In  the  IVP  for  the  integral  form 
of  the  conservation  laws,  initial  values  of  cell  averages  of  the  dependent  variables  will 
be  defined.  Given  such  initial  values,  the  three  steps  used  to  set  up  the  discretization 
procedure  for  the  integral  form  of  the  conservation  laws  Eire  as  follows. 

(i)  Define  dependent  variable  polynomials  in 

each  cell  so  that  the  cell  average  of  the  polynomial  approximation  matches  the  cell 
average  of  the  dependent  variable  which  is  either  given  as  part  of  the  initial  value 
specification  or  obtained  by  updating  the  cell  averages  during  subsequent  steps  of  the 
solution  process.  This  process  of  defining  pointwise  polynomial  behavior  from  known 
values  of  cell  averages  is  called  the  “reconstruction  procedure” .  At  quadrature  points 
on  boundary  faces  and  at  other  locations,  it  may  also  be  required  to  construct  the 
polynomials  by  fitting  them  to  known  pointwise  values  in  addition  to  cell  average 
values  of  the  dependent  variables. 

(ii)  Evaluate  the  polynomials  at  all  quadrature  points. 

This  will  lead  to  “left”  and  “right”  values  at  each  quadrature  point  which  lies  on  a 
face  common  to  two  cells. 

(iii)  Construct  the  solution  of  a  local  Riemann  problem  using  these  left  and  right  states 
and  from  this  evaluate  the  numerical  flux.  Use  such  numerical  fluxes  in  Eq.  26a  or 
Eq.  26b. 

(iv)  Steps  (i)-(iii)  complete  the  discretization  of  the  right  hand  side  of  Eq.  19.  To  the 
resulting  semi-discrete  system  of  equations,  we  apply  a  suitable  time-integration  pro¬ 
cedure. 

Features  of  Structured  Grid-Based  CEM 

•  Maxwell’s  equations  in  differential  conservation  form 

•  Explicit  Lax-Wendroff  upwind  finite-volume  scheme 

•  Multizone  structured  gridding 

•  Convenient  bookkeeping/data  structure 

•  Efficiently  vectorizable 

Disadvantages 

•  Multizone  gridding  process  is  tedious  and  time  consuming 

•  Zonal  interface  boundary  conditions  degrade  execution  efficiency 

•  Difficult  to  grid  certain  interior  regions 
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Featvires  of  Unstnictured  Grid-Based  CEM 

•  Maxwell’s  equations  in  integral  conservation  form 

•  Any  region  can  be  gridded 

•  Higher-order  basis  functions  for  solution  variables  in  each  cell 

•  Well  suited  for  MIMD  architectures 

•  Can  handle  special  regions  (thin  wire,  crack,  •  •  •)  through  choice  of  basis  functions 
Disadvantages 

•  Accuracy  may  be  an  issue  for  arbitrarily  arranged  unstructured  cells 

•  Vectorization  may  be  difficult 

Unstructured  Grid-Based  CEM  Issues 

•  Geometry /unstructured  grid  generation 

—  surface  definition 

—  surface  grid 

—  volume  grid 

•  Integral  equation  based  finite— volume  scheme 

—  basis  function  in  each  cell 

—  flux  evaluation  at  interface 

—  time  discretization 

•  Graphics/visualization 

In  the  development  of  an  unstructured  grid-based  CEM  method,  the  present  work 
addresses  some  of  these  issues. 
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4.0  CONSERVATION  FORM  FOR  MAXWELL^S  EQUATIONS 

We  can  regard  equations  (1)  and  (2)  as  a  first-order  system,  of  partial  differential 
equations  for  the  time  development  of  the  dependent  variables  B  and  D: 


dt\D)  \-V  xH )  \-j) 


(27) 


At  any  instant  the  functions  E,  H,  and  J  can  be  derived  from  the  current  values  of  D  and 
B,  and  by  substitution  into  Eq.  (27)  they  determine  new  values  for  the  time  derivatives. 

In  order  to  apply  CFD-based  conservation-law  form  finite-volume  methods,  Eq.  (27) 
is  rewritten  in  the  form  of  Eq.  (15),  where: 


Q=< 


0  ) 

(  D,/e 

-Dy/e 

0  ^ 

By 

-DJe 

0 

Dx/e 

0 

\b. 

D, 

^  ;  S  =  < 

Dyle 

0 

>  •  F  =  < 

1  1 

^  \G  =  < 

0 

By!  P 

^  ;  <S  =  < 

0 

-Jt 

Dy 

Bz/f^ 

0 

-Bx/p 

-Jy 

1  D,  J 

-By/p  J 

1  Bx/p  > 

\  0  > 

<.  z  ' 

(28) 

If  we  define  unit  vectors  x,y,z  along  three  orthogonal  coordinate  directions  in  space, 
equation  (28)  can  be  written  in  a  compact  form  as: 


xxE  yxE  \  df  zxE  0  \ 

dt\Dj^  dx\-xxHj'^  dy\-yxH  dz\-zxH)  \-J  J  ’ 


(29) 


in  which  every  one  of  the  six  component  equations  is  analogous  to  the  conservation  law 
dp/dt  +  V  •  J  =  0  which  holds  for  electric  charge.  In  a  formal  sense,  one  can  consider  the 
set  of  six  unknowns  (B,  D)  as  a  “vector”  density  Q,  and  define  a  “tensor”  flux  F  with 
components: 


F.  = 


X  X  E^ 
-X  X  H 


F  =  (  ^  ^ 

"  [-yxH 


F.  = 


z  X  E 
-zxH 


(30) 


such  that  equation  (29)  becomes: 

dQ/dt  +  V  •  F  = 


0 

-J 


^5 


(31) 


This  rather  abstract  representation  as  a  “system  of  conservation  laws”  provides  a 
convenient  starting  point  for  the  development  of  integration  algorithms. 

Also,  one  can  transform  the  equations  to  arbitrary  curvilinear  coordinates 
without  alteration  to  this  basic  form: 


^  \  j)  V-f  X  H/jJ 


drj 


T]  X 


E/J  \ 


-ifxH/j)  di:\-CxH/jJ  J 


EU  \_S_ 


(32) 
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where  J  is  the  Jacobiaa  of  the  coordinate  transformation  and,  for  example, 

In  addition  to  the  differential  form,  one  can  use  the  divergence  theorem  on  any  volume 
V  to  obtain  the  integral  form  of  MaxwelU’s  equations: 


—n  X  H  J 


■JU 


S  dV 


where  the  six  components  of  •  n  in  Eq.  (19)  are  (n  x  E,  n  x  H). 

Equation  (33)  replaces  all  the  spatial  derivatives  with  area  integrals  over  the  surface 
of  F. 

4.1  Interface  Flux  Forms 

In  order  to  evaluate  the  fluxes  f,  E,  and  Q  required  by  a  finite-volume  procedure, 
Eq.  (24),  at  appropriate  cell  boundaries,  first  the  characteristics  of  Eq.  (32)  are  an^yzed 
in  each  (t,0,  (^,»7),  and  (^>0  plane.  For  example,  the  characteristics  (eigenvalue^)  injhe 
(r,0  plane  are  given  by  solving  the  matrix  equation  1-4  -  XI\  =  0,  where  A  =  dS/dQ  is 
the  Jacobian  matrix. 

For  a  locally  isotropic  medium,  D  =  ^  and  B  =  ,  the  matrix  A  has  six  eigenvalues, 

two  vanishing  (Ai,2  =  0),  two  equal  tQj:||‘l,  and  two  equal  to  -c\(  |,  where  c  =  1/^  and 
l^j  =  /^2  Referring  to  Fig.  4.1.1,  the  negative  and  positive  eigenvalues  are 

indicated  at  a  cell  boundary  (interface).  The  intent  is  to  compute  the  interface  flux  given 
the  left  state  —  and  the  right  state  +  on^ither  jide  of  the  interface.  For  a  linear 
system,  across  a  characteristic  A,  the  variation  in  Q  and  £  are  related  by  a  jump  condition 
-  Ale  I  +  |F  1  =  0,  where  |  1  denotes  jump.  To  make  the  treatment  more  general,  a  resistive 

sheet  is  introduced  at  a  cell  boundary  (Fig.  4.1.1)  whose  resistivity  is  denoted  by  O-lcrd) 
where  a  is  the  conductivity  and  d  is  the  thickness  of  the  sheet.  The  resistive  sheet  can 

allow  for  total  tangential  magnetic  fields  (j  x  to  jump  across  the  sheet.  Corresponding 
to  the  ^-direction  interface,  the  following  jump  relationships  are  written^^. 

Across  —  Characteristics  ^A  =  —c\^\ ,  ^ 

{B*  -  B-)  = X  {E*  -  E-) 


Across  +  Characteristics 

( 5+  -  B”)  =  ^  X  (E+  -  E**) 

At  the  Interface  (A  =  0)  (Allowing  for  a  Resistive  Card) 
fx(E**-E*)  =  0 
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f  X  {H**  -  ^* )  =  -(Td  X  f  X  E*)  (34d) 

Relationship  Eq.  (34d)  is  obtained  by  integrating  Eq.  (32)  over  an  infinitesimal  strip  area 
enclosed  by  a  contour  containing  the  resistive  sheet.  Solving  the  above  four  relationships, 
one  can  write  the  expressions  for  the  interface  fluxes: 


f  X  E*  =  f  X 


(xH*  =(x 


(xH** 


{ [f)+(€c)+  +  f  X  ff+]  +  [(ec)-B-  -  f  X  if-] } 

(ec)“  +  (ec)+  +  crd 

I  [l  +  ad{fic)+]  {H-(fic)-  +  f  X  E-)  +  [(pc)+if+  -  f  x  E+]  | 
(/ic)+  +  (pc)-  +  o-d(pc)+(pc)- 

{ [1  +  (Td(pc)-]  {(fic)+H+  -  f  X  E+)  +  [(pc)-E-  +  f  X  E-] } 
(pc)+  +  (pc)-  +  ad(pc)-(pc)+ 


The  superscripts  +  and  —  denote  the  right  and  left  state  at  a  cell  interface.  When  ad  =  0, 
the  total  tangential  electric  and  magnetic  fields  axe  continuous  across  a  material  interface. 
When  ad  -+  oo,  the  resistive  sheet  represents  a  perfectly  conducting  surface,  and  the 
total  tangential  electric  field  (  x  E*  goes  to  zero  satisfying  Maxwell’s  boundary  conditions 
for  a  perfectly  conducting  surface  (to  be  discussed  later).  The  corresponding  tangential 

magnetic  fields  x  depend  only  on  E~  and  H~,  and  similarly  x  depends 

only  on  E"*"  and  E"*".  Equation  (35)  allows  for  the  material  properties  (e,  p)  to  jump  any 
amount  at  an  interface.  In  free  space  where  c"^,  c”,  p"^,  p-,  and  e~  are  normalized  to 
unity,  Eq.  (35)  reduces  to  a  simpler  form. 


Given  a  left  (— )  state  and  a  right  (+)  state,  the  process  of  computing  the  intermediate 
states  (*  and  **)  is  usually  referred  to  as  the  “Riemann”  problem. 


Additional  information  on  these  interface  flux  forms,  and  on  the 
eigenvalue/eigenvector  structure  of  the  Maxwell  equations  is  given  in  Section  4. 

4.2  Characteristics  of  Maxwell’s  Equations 


The  characteristics  of  Eq.  (21)  are  analyzed  by  considering  the  Jacobians 
1^,  and  For  example,  the  Jacobian  matrix  ^  for  Maxwell’s  equations,  Eq.  (32), 
becomes: 

0 


A-^ 

dQ 


0 

0 

0 

0 

-^z/e 

iyh 


0 

0 

0 

0 

-izh 


0 

0 

ixh 

0 


0 

0 

0 

0 


-iz|^J■ 

0 

^z!  p 
0 
0 
0 


iy!  M 
iziti 
0 
0 
0 

0  J 


(36) 


which  can  be  symbolically  written  as: 


BQ  V(fA)x  0  ) 


(37) 
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|x{E**-E*)  =  0 

=(odlE* 


Fig.  4.1.1.  Interface  Flux  Treatment  with  a  Resistive  Sheet. 


where: 


0 

-iziii 

0 

Cxifi 

V^/ 

1 — 

H 

1 

0 

(38) 


The  eigenvalues  of  the  Jacobian 


matrix  A  in  Eq.  (36) 


are  obtained  from  solving: 


|yl-AJ|  =  0 


(39) 


and  the  right  eigenvector  corresponding  to  a  A:th  eigenvalue,  Ajt,  is  obtained  by  solving: 


\A-Xkl\{rk]  =  0 


(40) 


For  Maxwell’s  equations,  solution  to  Eq.  (39)  reuslts  in  six  eigenvalues:  two  vanishing, 
two  equal  to  c|if|,  and  two  equal  to  — c|^|,  where  c  =  l/y/jle,  ^  =  (G,^j,i^z),  and  = 

Defining  the  upper  and  lower  halves  of  the  eigenvector  f  as  e  and  h,  respectively,  so 
that  f  =  one  finds  corresponding  to  the  zero  eigenvalues: 


(41) 


where  ^  =  ^/||].  Corresponding  to  the  positive  or  negative  eigenvalues,  one  finds  that 
e,  h,  and  ^  form  a  mutually  orthogonal  system  with  h  =  ±(^  x  e/ec),  and  where  the 
upper  sign  goes  with  the  negative  eigenvalue.  For  example,  corresponding  to  the  two 
negative  eigenvalues  one  can  construct  the  two  independent  set  of  eigenvectors  by  choosing 


19 


c  =  a  tangent  vector  in  the  direction  »/  on  a  constant-^  plane  for  one,  and 

e  =  a  tangent  vector  in  the  direction  ^  for  the  other.  Once  e  is  chosen  h  can 

be  obtained  using  the  mutually  orthogonal  property  h  =  x  e/ec). 

For  hyperbolic  systems  of  equations,  the  Jacobian  has  reed  eigenvalues  and  a  linearly 
independent  set  of  eigenvectors^^’^^.  The  equations  are  also  characterized  as  “linear”  if: 

VcAfcTjt=0  (42) 

and  “nonlinear”  if: 

VfiA*  •  Tik  ^  0  .  (43) 

It  can  be  verified  that  time-domain  Maxwell’s  equations  are  hyperbolic  and  linear.  For  a 
linear,  hyperbolic  equation  of  the  form 

9t  +  /i  =  0  (44) 


one  can  write  the  jump  condition^"*: 


-A[q]  +  [/]  =  0  (45) 

where  [q]  and  [f]  are  the  changes  across  an  eigenvalue  A.  This  jump  condition  is  used  in 
deriving  the  interface  flux  forms  in  Section  3.2. 

For  a  nonlinear  system,  the  jump  condition  is  defined  only  across  a  surface  of  discon¬ 
tinuity  S{x,t)  =  0,  such  that: 

St[q]  +  SM]  =  0  (46) 

where  (St/Sx)  =  —{dx{dt)  =  A  is  the  characteristic  speed  of  the  discontinuity.  For  a  lineeir 
system  A  is  replaced  by  A,  the  eigenvalue. 


4.3  The  Riemann  Solver;  Preliminaries 


The  wave  solutions  of  Maxwell’s  equations  can  be  derived  in  a  different  way  by  looking 
for  characteristic  combinations  of  E  with  H  (or  D  with  B)  that  propagate  without  change 
along  a  particular  direction  in  space.  The  form  of  these  combinations  can  be  illustrated 
using  one  space  dimension  (say,  i),  in  which  case  Maxwell’s  equations  reduce  to: 


dt  dx  '  dt  dx 


(47) 


In  free  space,  Ey 
matrix  form  as: 


or: 


Dy/eo  and  Hz  =  Bz/fio,  so  the  two  equations  can  be  rewritten  in 


dt 


Bz 


+ 


dx 


(°o-'“o)(d‘)-°  ’ 


■^Q  +  =  0,  where  Q  =  {Bz,Dy) 


(48a) 

(486) 
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By  direct  substitution  into  (47),  it  is  easy  to  show  that: 


where  c  =  lly/jl^  is  the  velocity  of  light  in  vacuum. 

These  are  just  two  scalar,  first-order  equations  having  the  general  solutions: 


Bz  +  —  gi(^x  —  ct)  , 

(50a) 

V 

Bz- ,[^Dy  =  g2{x  +  ct)  , 

V  ^0 

(506) 

the  first  representing  a  wave  traveling  toward  increasing  values  of  x  and  the  second,  a  wave 
traveling  in  the  opposite  direction.  The  quantity  =  \/ is  known  as  the  free— space 
impedance. 

From  (50)  the  most  general  solution  of  the  one-dimensional  Maxwell’s  equations  can 
be  written  as: 

B,=^gi{x-ct)  +  g2ix  +  ct)  (51a) 

Dy  =  (a:"-  ct)  -  g2{x  +  ct)]  /Co  •  (51^) 

So,  whatever  initial  distribution  of  B  and  D  we  are  given  at  <  =  0,  we  can  decompose  it 
uniquely  into  a  right— moving  and  a  left— moving  wave,  each  of  which  travels  at  the  velocity 
of  light.  In  the  (i,t)  plane,  the  right-moving  combination  gi  =  B  +  Co-D  is  seen  to  be 
constant  along  the  characteristic  lines  x  —  ct  =  const.,  while  g2  =  B  —  CoD  is  constant 
along  the  lines  x  +  ct  =  const.  Thus,  the  value  of  B  at  any  point  in  the  upper  {x,t)  plane 
can  be  constructed  as  shown  in  the  figure  from  the  initial  values  of  B  and  D  at  xq  =  x  —  ct 
and  x\  =  X  +  ct. 


Fig.  4.2.1  Right  and  left  moving  characteristics. 
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Another  way  to  view  the  characteristic  combinations  gi  =  B  +  ^qD  and  g2  =  B~(oD 
is  as  defining  eigenvectors  of  the  2x2  matrix  A  in  equation  (48): 


(B,D)  =  (1,+Co  for  right  moving  waves  (52a) 


(jB,  D)  =  (1,  —Co  for  left  moving  waves 

U'  (-i-) 


(526) 

(53a) 

(536) 


This  reveals  a  fundamental  connection  between  the  abstract  representation  as  a  system  of 
conservation  laws  and  the  propagation  characteristics  of  the  original  problem:  the  eigenvec¬ 
tors  of  the  abstract  system  (known  as  the  Riemann  invariants)  provide  a  way  to  decompose 
the  initial  data  that  allows  us  to  construct  the  solution  at  all  later  times. 


For  a  linear  hyperbolic  equation  of  the  form  qt  +  fx  =  0,  one  can  diagonalize  the 
Jacobian  matrix  A  =  ^  such  that  A  =  where  L  is  the  left  eigenvector  matrix, 

R  is  the  right  eigenvector  matrix,  and  [A]  is  the  diagonal  matrix  made  up  of  eigenvalues. 
Multiplying  throughout  by  L,  one  can  write  the  equation  {Lq)t  -|-  [A](L?)i  =  0,  where  now 
the  combination  (Lq)  is  the  Riemann  invariant. 

4.4  Discretization  and  Dissipation 

The  true  initicd-value  problem  we  must  solve  numerically  is  to  determine  the  behavior 
of  the  electromagnetic  fields  on  a  discrete  mesh,  from  data  originally  prescribed  on  this 
same  mesh.  In  one  dimension,  the  fields  B{x,t)  and  D(x,t)  are  represented  by  values 
6,(t)  and  d,(t)  attributed  to  each  interval  of  the  space  coordinate  x  that  lies  within  the 
computational  domain  (say,  0  <  x  <  L).  In  our  finite-volume  approach,  these  values  can 
be  thought  of  either  as  approximating  the  true  values  of  B  and  D  at  the  center  of  each 
interval  at  a  given  time,  or  as  the  averages  of  B  and  D  over  the  interval.  More  generally, 
B  may  be  evaluated  on  a  different  spatial  mesh  from  D,  or  the  behavior  of  B  and  D  within 
each  interval  may  be  modeled  by  high-order  polynomials. 

In  any  case,  one  cannot  directly  solve  the  partial  differential  equations  for  B  and  D, 
but  instead  must  seek  a  discrete  set  of  equations  for  {6,}  and  {d,}  whose  solution  for  any 
later  time  approaches  that  of  the  original  Maxwell’s  equations  as  the  mesh  is  refined.  This 
is  the  basic  concern  of  all  numerical  integration  techniques.  In  what  follows,  we  describe 
a  particular  approach  to  this  problem  that  uses  the  Riemann  invariants  to  achieve  both 
accuracy  and  stability  in  the  integration  process. 

This  method  was  developed  primarily  by  S.  Osher^^  and  various  co— workers,  includ¬ 
ing  Sukumar  Chakravarthy  of  the  Rockwell  Science  Center.  They  considered  the  general 
hyperbolic  system: 

+  =  0  .  (54) 
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where  w{x,t)  is  a  vector  of  unknown  fields  and  the  fiux  f{w)  is  a  vector-valued  function 
whose  Jacobian  matrix  df/dw  has  only  real  eigenvalues. 

Discretizing  the  space  variable  in  this  equation  forces  one  to  replace  the  spatial  deriva¬ 
tive  with  an  approximation  based  on  finite  differences.  The  order  of  this  approximation, 
in  the  sense  of  Taylor  series  expansion,  is  then  found  to  determine  the  order  of  accuracy 
of  the  computed  fields.  Furthermore,  the  details  of  this  approximation  place  limits  on 
the  stability  of  any  scheme  chosen  to  integrate  (54)  forward  in  time.  In  particular,  we  do 
not  want  the  ragged  form  of  the  discretized  initial  data  u;,(0)  to  give  rise  to  increasing 
raggedness  in  the  solution  tuj(t)  8>t  later  times. 

What  Osher  did,  in  essence,  was  to  rewrite  (54)  in  terms  of  its  Riemann  invariant 
combinations,  and  determine  for  each  invariant  an  appropriate  discretized  flux.  Just  as  in 
our  ID  Maxwell  example,  this  flux  carries  information  only  from  the  past  locations  of  the 
wave  to  the  present,  md  it  is  therefore  called  an  “upwind  approximation. 

The  first  step  in  Osher’s  procedure  is  to  find  the  eigenvectors  r*  of  dfidw,  together 
with  their  associated  eigenvalues  At.  As  we  have  seen  in  the  example,  these  eigenvectors  are 
the  combinations  of  the  components  of  w  that  propagate  locally  as  simple  waves  gk{x—\kt)‘‘ 


dw 


Tk  =  Afcr/t 


(55) 


For  the  special  case  of  a  linear  dependence  of  /  on  u;,  we  note  that  the  matrix  dffdw  is 
constant,  and  the  eigenvectors  and  ei^n values  are  independent  of  w. 

The  second  step  is  to  construct  a  two-point  approximation  h+{wj,Wj+i)  for  the  flux 
entering  the  jth  cell  from  the  interface  between  cells  j  and  ;  -f  1.  Osher  does  this  by 
finding  a  value  that  is  consistent  with  data  propagated  from  the  interior  of  both 

neighboring  cells  ana  putting: 

/i+(ujj,u;j+i)  =  /(u;*+i/2)  •  (56) 

A  similar  construction  at  the  [j  —  l,j)  interface  yields  a  value  consistent  with  data 

in  cells  ;  —  1  and  j  and  an  approximation: 

for  the  flux  entering  the  jth  cell  from  the  j  -  1/2  interface.  These  approximations  can  be 
used  directly  to  obtain  a  first-order  accurate  approximation  to  (54): 


dwi  1 

"X  +  A 


/(“'I+1/2)  - 


=  0 


(58) 


More  importantly,  the  first-order  fluxes  form  the  basis  for  constructing  integration  schemes 
of  higher  order,  such  as  the  Lax-Wendroff  procedure  used  in  this  report. 

How,  then,  are  these  interface  values  of  w  constructed  from  the  Riemann  invariants? 
In  Osher’s  method  an  integration  is  performed  in  the  abstract  space  of  the  possible  values  of 
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w,  starting  (in  the  first  case)  from  Wj  and  ending  at  tWj+i,  and  following  a  path  determined 
by  those  eigenvectors  that  represent  waves  propagating  toward  from  the  interface  from 
either  side.  The  first  leg  of  this  integration  puts: 


dw 

dsi 


ri(u)(si)) 


(59) 


where  t\  is  the  eigenvector  that  corresponds  to  the  wave  traveling  at  the  greatest  speed 
Ai  <  0^  away  from  the  interface  from  the  interior  of  the  jth  cell.  The  first  leg  ends  at 

a  point  Wj^i  which  we  have  to  leave  imdetermined  until  the  path  finally  reaches  wj^i.  The 
second  leg  puts: 

^  =  r2{w{s2))  (60) 

and  so  forth,  until  all  the  waves  propagating  from  j  +  1/2  to  j  have  been  used.  The  (as 
yet  undetermined)  value  of  w  at  this  end  point  is  Wj^ij2‘  Now  the  waves  for  which  A*;  =  0 
axe  used  in  the  same  fashion,  i.e.,  one  integrates 


—  =  rk{w{sk))  (61) 

along  each  such  eigenvector.  The  end  point  of  this  series  of  integrations  will  be 

The  last  series  of  legs  involve  the  waves  that  propagate  toward  the  interior  of  cell  j'  + 1 
from  the  j  +  1/2  interface,  done  in  the  order  from  smallest  speed  to  largest.  Finally,  on 
the  last  leg  one  has: 

^  =  r,rx{w{s„,))  ,  (62) 

dSjjx 

and  the  final  end  point  is 

The  m  known  components  of  Wj  and  wj+i  provide  just  enough  information  to  deter¬ 
mine  the  m  scalar  quantities  {s*}  eind  therefore  the  interface  values 

As  an  important  illustration  of  this  procedure,  consider  the  one-dimensional  Maxwell’s 
equations  that  describe  an  inhomogeneous  medium,  that  is,  a  medium  in  which  e  and  fi 
are  functions  of  x: 


If  we  take  w  to  be  simply  {B,  D)  then  this  equation  is  not  in  the  form  of  (54)  because  the 
flux  depends  explicitly  on  x  through  e  and  /x.  To  get  around  this  difficulty,  we  define 
and  as  additional  dependent  variables  e  and  m,  and  write: 


dt 


=  0 


(64) 
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Then  the  Jacobian  matrix  becomes: 

(0  c  D  0 

m  0  0  B 

0  0  0  0 
0  0  0  0 

■  which  has  two  zero  eigenvalues  A2  =  A3  =  0,  one  positive  eigenvalue  A4  =  c  =  em 
negative,  Ai  =  — c.  The  corresponding  eigenvectors  axe: 


ri  = 


rj  = 


rs  = 


r4 


Along  the  first  leg  of  integration,  we  now  have: 


dB  dD  . 

ds  ~  ^  ^  ds  ^  '  ds  ds 

■  Bj  —  cs  ;  ~  ~  ~  ^ 


B 


j+1/2 


ruj 


■®jVl/2 


Along  the  second  leg: 


dB  .  dD  ^  de  dm 

17  =  “  '  17  =  ^  ’  57  =  -'  ’  17  =  “ 


Bj,2  —  Bj^^i2  ;  In 


■^i+l/2 


=  s  =  —  In 


fill 


mj,2  =  nxj 


^j,2-Dj,2  —  ^j^j+1/2  "*■  ^j+1/2  ~  ^j+i/2  • 


Along  the  third  leg: 


dB  „  dD  .  de  .  dm 

17  =  -®  ’■  17  =  “  ■  17  =  “  '  17  =  -'" 


ft-  — m-o  ’  ~i+l/2 

-Dj,2  nij^2 

mj,zB*Xii2  =  ^i^jVi/2  -^>+1/2  =  ■^jVi/2 


and  along  the  final  leg: 


ds 


=  e 


ds 


de 


dm 


’  !  ds  ^  ds 


;  Cj+j  —  ej^  =  0  =  mj+i  —  mj,3 


(66) 

and  one 

(66) 

(67) 

(68) 

(69) 

(70) 

(71) 

(72) 

(73) 

(74) 

(75) 

(76) 


25 


(77) 


(Si+i  -  =  i  -  oy,,,) 

Putting  these  relations  together,  one  finds  finally: 


j+i/2  J  j+i/2  (cj/mj)  +  (cj+i/mj+i) 


(78a) 


■^i+l/2  -  ^j^j+1/2 


_  “t"  "1“  [^i+i-^i+1 

(Cj/Cj)  +  (Cj+i/Cj+i) 


(786) 


We  thus  have  found  that  the  appropriate  fluxes  involve  a  weighted  average  of  the  Riemann 
invariants  on  either  side  of  the  interface.  The  electric  terms  that  now  appear  in  H* ,  and 
the  magnetic  terms  in  E*,  when  these  expressions  are  substituted  into  (58),  result  in  a 
diffusion-like  second  difference  term  that  acts  to  damp  the  integration  process.  These  terms 
in  fact  provide  the  necessary  stability  that  is  missing  from  the  simple  central  difference 
approximation. 


The  generalization  of  this  approach  to  three  dimensions  is  straightforward:  The  flux 
contribution  along  each  of  the  three  coordinate  directions  is  computed  independently,  and 
then  all  three  flux  terms  Eire  added  to  give  dwj dt  in  first  order.  Along  the  ^  direction,  the 
fluxes  are: 


f  X  E*  =  f  X 


[c-p-  -  I  X  H-]  +  [C+.D+  + 1  X  H+] 
(ce)-  +  (ce)+ 


(79a) 


f  X  f  X 


[c-B-  +  I  X  E-]  +  [c+i+  -  X  E+] 
(fic)-  +  (/ic)+ 


(796) 


which  are  the  same  as  the  fluxes  given  by  equation  (35)  when  there  is  no  resistive  sheet  at 
the  cell  interface. 
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S.O  UNSTRUCTURED  FINITE-VOLUME  TREATMENT 
5.1  Space/Time  Discretization 

The  major  feature  of  our  discretization  approach  that  distingtiishes  it  from  other 
finite~volume  and  finite— difference  procedures  is  that  the  electric  and  magnetic  field  un¬ 
knowns  are  co-located  in  both  space  and  time,  rather  than  being  assigned  to  two  inter¬ 
penetrating  spatial  grids  and  separated  a  half— step  in  time.  These  field  unknowns  are  the 
volume  averages  of  E  and  H  within  each  cell  in  the  space-filling  grid. 

Staggered-grid  methods  automatically  achieve  second-order  accuracy  in  space,  while 
co-located  field  algorithms  require  near-neighbor  corrections.  However,  there  is  a  funda¬ 
mental  equivalence  between  these  methods  in  terms  of  the  achievable  accuracy  and  stability 
of  the  integration  process. 

Both  approaches  typically  use  explicit  time  integration,  which  means  that  the  upper 
limit  on  the  allowable  size  of  the  time  step  At  is  determined  by  the  physical  size  and  shape 
of  the  smallest  cells,  corresponding  roughly  to  the  time  that  light  takes  to  cross  one  of 
these  cells.  Implicit  integration  schemes  can  choose  larger  time  steps,  but  they  require  the 
inversion  of  a  banded  matrix  the  size  of  the  whole  grid,  and  their  ability  to  preserve  phase 
information  is  not  known. 

The  unstructured  algorithm  we  have  developed  is  applicable  to  any  grid  that  fills 
the  computational  domain  with  polyhedral  cells.  In  particular,  we  have  implemented  the 
necessary  bookkeeping  procedures  to”  deal  with  hexahedra  (such  as  cubes),  tetrahedra, 
and  prisms  (translations  of  a  triangle  out  of  its  plane  of  definition).  A  generalization  to 
curved  interfaces  is  straightforward,  but  not  regarded  as  useful  at  this  time.  These  three 
polyhedral  geometries  allow  us  to  specialize  the  code  in  various  ways,  for  instance  to  check 
against  the  structured-grid  Lax-Wendroff  algorithm  in  our  present  solver  RCS3D,  or  to 
run  purely  two-dimensional  cases. 

Each  polyhedron  in  the  computational  domain  is  specified  by  the  location  (i,j/,z)  of 
its  vertices  in  physical  space.  From  these  locations,  all  the  necessary  geometrical  quantities, 
including  areas,  surface  normals,  and  centroidal  locations  are  computed.  As  stated  earlier, 
each  field  unknown  attributed  to  a  given  polyhedron  is  considered  to  be  an  averap  of  the 
field  over  the  volume  of  the  polyhedron.  The  six  components  of  E  and  H  at  one  time  level 
are  thus  stored  according  to  an  index  a  that  runs  over  all  polyhedra.  Quantities  related 
to  the  polyhedral  faces,  such  as  face  normals,  are  stored  according  to  another  index  that 
runs  over  all  faces. 

The  interior  faces  of  a  given  polyhedron  are  kept  distinct  from  those  of  the  neighboring 
polyhedra  that  share  faces  with  it  in  a  purely  geometrical  sense.  This  allows,  for  instance, 
for  any  type  of  impedance  boundary  condition  to  be  applied  at  the  boundary  between 
cells.  Thus,  each  polyhedral  face  has  a  co— face  with  a  distinct  face  index,  and  each  such 
co-face  belongs  to  its  own  polyhedron. 
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5.2  Polynomial  Representation  and  Least  Squares 

To  go  beyond  representation  of  the  fields  as  simple  voliime  averages,  we  have  chosen 
initially  to  implement  linear  polynomial  functions  for  both  E  and  H,  Higher  order  poly¬ 
nomial  representations  will  follow  the  same  general  procedures.  The  essential  question  is 
how  the  higher  order  terms  in  these  polynomials  axe  to  be  determined  from  near-neighbor 
data,  so  as  to  achieve  the  desired  level  of  accuracy  within  each  cell.  In  our  imstructured 
approach,  this  evaluation  is  closely  tied  to  the  time  integration  procedure  through  the  Rie- 
mann  fluxes  at  each  interface.  Ultimately,  this  preserves  time  accuracy  as  well  as  accuracy 
in  space. 

In  the  first  step  of  this  method,  first-order  Riemann  fluxes  are  constructed  at  each 
interface  of  a  cell  from  the  volume-averaged  fields  on  either  side  of  the  interface.  For 
Maxwell’s  equations,  these  fluxes  are  the  tangential  field  components  just  inside  the  bound¬ 
aries  of  the  cell.  To  complete  the  specification  at  the  cell  surface,  the  normal  components 
of  E  and  H  are  taken  to  be  the  normal  components  of  the  volume-averaged  fields.  This 
maintains  overall  charge  conservation  within  the  cell. 

These  boundary  data  are  sufficient  to  determine  all  the  terms  in  a  linear  polynomial 
fit  to  either  E  or  H  hy  &  procedure  such  as  least-squares  minimization  of  the  fitting  error 
integrated  over  the  boundary.  If  we  denote  the  vector  polynomial  to  be  fitted  as  A  and  its 
boundary  values  as  A* ,  then  the  quantity  to  be  minimized  is 


/, 


(cell  boundary) 


dS 


(80) 


Taking  derivatives  of  e  with  respect  to  each  polynomial  coefficient  in  A  results  in  a  non¬ 
singular  set  of  linear  equations  for  these  coefficients.  For  consistency,  one  constrains  the 
constant  term  in  A  to  be  equal  to  the  known  volume  average  of  A  over  the  cell. 

A  separate  set  of  equations  is  obtained  for  each  Cartesian  component  of  A.  If  one 
writes,  e.g., 


A  ,  /  X  dAx  ,  X  ,  .  dAx 

Ax{r)=  {Ax)a  +  (a:  -  Xa)  +  (y  -  y^)  +  {z  -  z^) 


(81) 


where  the  angular  brackets  denote  volume  averaging  and  e.g.,  Xa  =  (x)q,  then  these 
equations  become 


(x-Xa)  {x  -  Xca)^^  +  {y  -  +  (z  -  Za) 


dx  '  dy 

(x  -  Xa)[Al  -  {Ax)o,]dS 


do 

Jq,  (2/  -  2/“)  [(^  "  ^  +  (2/  -  Vo)  ^  +  (^  - 


dAx 

dz 


dAi 


dS 


(82a) 


dS 


=  [  iy-yo)[Al-{Ax)a] 

Jda 


(826) 


dS 
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(82c) 


=  f  (z-  Za)  [A*  -  (A*)a]  dS 
Jda 


sdA, 

^  dz 


dS 


> 


where  we  have  denoted  the  cell  boundary  as  doc.  These  equations  can  be  solved  by  inverting 
the  matrix  M  whose  elements  are  the  quadratic  moments 


Mij  = 


ra)  j  •  {r-ra)dS 


(83) 


where  t  and  j  are  unit  vectors  in  the  respective  coordinate  directions. 

For  a  linear  polynomial  fit,  there  is  a  simpler  alternative  procedure  that  we  have 
implemented  to  evaluate  these  linear  terms.  From  the  divergence  theorem,  the  average 
value  of  any  derivative  over  the  cell  voltune  can  be  rewritten  as  a  surface  integral; 


(84) 


where  h  is  the  unit  outward  normal  on  the  boundary  da  and  Va  is  the  cell  volume.  In 
particular,  if  p  is  a  linear  function  of  r  then  dpfdx  is  constant  and  equal  to  this  volume 
average,  which  can  be  calculated  just  from  the  values  of  p  on  the  boundary.  For  every 
component  of  A,  we  can  replace  its  l^undary  values  by  the  corresponding  component  of 
A*  to  obtain  the  approximation 

f  nA*dS  =  ,  (85) 

Va  Jda 


which  is  equivalent  to  using  n  as  the  weight  in  the  method  of  weighted  residuals  applie_d  to 
the  difference  A -A*.  The  quantity  Ka  is  a  vector  dyadic.  Since  we  have  chosen  h- A*  = 
h  •  {A)a^  we  can  make  use  of  the  vector  identity  a  =  n  (n  •  a)  —  n  x  (n  x  a)  to  rewrite  the 


integral  as 


Ka  =  f  nnxjnx  ((A)a  -  A*)} 
Jda  ’■ 


dS 


(86) 


which  will  be  more  convenient  to  compute  in  terms  of  the  tangential  components  of  A*. 
This  particular  weighting  can  be  shown  to  result  from  a  variational  principal  that  assumes 
each  Cartesian  component  of  A*  is  the  boundary  value  of  a  solution  of  Laplace’s  equation 
inside  the  cell. 


5.3  The  Unstructured  Second— Order  Algorithm 

An  algorithm  that  maintains  second-order  accuracy  in  both  space  and  time  can  be 
constructed  from  the  linear  polynomial  representation  as  follows: 

-  #•  /  n-F(QT)dS  (87) 

ZVa  Jda 
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Jf?  =  i  /  hQTiS^^f  ft  (n  X  {n  X  [W):  -  Or] })  iS  (88) 

Vet  Jda  Jda 

Q^+'l^(f)  =  {Q)"*'/^-l(r-r„)-K"  forfinceUa  (89) 

(«)?^‘  =  (0)r-^  /  •  (90) 

Here  we  have  written  Maxwell’s  equations  symbolically  as 

^+V-F(Q)  =  0  ,  CI=[B,b)  ,  (91) 


and  the  solution  of  the  Riemann  problem  just  inside  a  cell  interface  is  denoted  Q*.  This 
solution  depends  only  on  the  values  of  Q(f)  immediately  on  either  side  of  the  interface. 
These  are  the  cell-average  values  for  Q*"*  and  the  hnear  polynomial  values  for 
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6.  UNSTRUCTURED  GRIPPING  AND  RESULTS 


There  are  three  steps  involved  in  the  unstructured  grid  setup.  1)  Gwmetry  modeling, 
2)  surface  grid  setup,  and  3)  volume  grid  setup.  Appendix  A3,  sections  of  a^  contract 
report  on  geometry  and  gridding  prepared  by  the  Computational  Fluid  Dynamics  group 
at  Rockwell  Science  Center,  include  a  detailed  writeup  on  the  unstructured  grid  generation 
process  using  a  procedure  called  the  “advancing  front  method.”  Once  the  unstructured 
grid  is  set  up  (hedahedron,  prism,  or  tetrahedron),  a  preprocessor  is  used  to  define  all  the 
quantities  required  by  the  unstructured  grid-based  Maxwewll’s  solver.  Some  of  them  are, 

1.  Cell  volumes  and  cell  face  normals, 

2.  Face/coface  connection  (from  this  information  one  constructs  the  face  and  node  near 

neighbor  connection),  and 

3.  Boundary  face  numbers  (both  body  and  outer  boimdary). 

This  grid  related  information  is  used  in  the  flux  summation  process. 

Problems  in  CEM  involve  arbitrarily  shaped  three-dimensional  geometries  that  need 
to  be  represented  properly  in  the  computer  simulation.  In  addition  to  the  external  shape, 
CEM  also  requires  modeling  the  interior  of  the  penetrable  structure.  Depending  on  the 
formulation  (differential  or  integral),  one  may  choose  either  a  structured  grid  or  an  un- 
structured  grid  setup. 

Two  gridding  issues  that  need  to  be  addressed  in  EM  computations  are;  1)  num¬ 
ber  of  grid  points  per  wavelength  to  properly  represent  the  fields  in  and  around  a  scat- 
terer;  and  2)  how  far  should  the  outer  boundary  be  placed  from  the  scattering  object 
to  adequately  simulate  the  nonreflecting  boundary  condition.  In  general,  the  number  of 
points/wavelength  is  not  determined  by  wavelength  alone,  and  involves  the  body  dimen¬ 
sions  (characteristic  body  size  with  respect  to  wavelength)  also.  The  outer  boundary 
location,  theoretically,  can  be  right  on  the  body  surface  itself;  however,  the  computational 
implementation  of  nonreflecting  boundary  conditions  requires  the  outer  boundary  at  a  few 
(2  to  5)  wavelengths  away  from  the  surface.  Again,  if  one  can  construct  higher  order  ac¬ 
curate  implementations  of  nonreflecting  boundary  conditions,  the  outer  boundary  c^  be 
brought  very  close  to  the  scattering  surface.  In  general,  the  necessary  grid  resolution  is 
provided  only  around  and  near  the  body  surface.  Between  the  body  and  the  outer  bound¬ 
ary,  the  mesh  is  allowed  to  stretch  resulting  in  very  crude  (3  to  5  points  per  wavelength) 
meshes  near  the  outer  boundary  regions. 

The  free  space  wavelength  is  reduced  to  smaller  values  inside  a  material  (as  e  and  /i 
become  large,  the  speed  of  propagation,  c  =  ^,  goes  down,  causing  the  wavelength  to 
scale  accordingly).  Thus,  the  grid  resolution  must  take  into  account  material  properties 
to  adequately  resolve  the  fields  inside  material  zones. 

The  number  of  grid  points  per  wavelength  required  depends  on  the  order  of  accuracy 
of  the  numerical  scheme.  A  second-order  accurate  scheme  usually  requires  at  least  ten  grid 
points  per  local  wavelength.  One  may  be  able  to  use  a  higher  order  scheme  and  minimize 
the  number  of  grid  points.  However,  as  the  order  of  accuracy  goes  up,  the  scheme  will  also 
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require  more  computations  per  grid  point,  which  may  offset  the  execution  savings  with 
fewer  grid  points. 

The  requirement  that  the  fields  are  resolved  accurately  with  proper  grid  resolution 
makes  CEM  problems  computationally  intensive,  requiring  large  scale  supercomputing. 
For  example,  to  compute  the  radar  cross  section  of  a  typical  aircraft  at  1  GHz,  even  if  one 
used  10  grid  cells  per  wavelength,  it  will  require  tens  of  millions  of  grid  points. 

Some  of  the  salient  features  of  the  current  CEM  capabilities  are 

•  Time-domain  Maxwell’s  equations 

•  Proven  algorithms  from  CFD 

•  Single  pulse  (multiple  frequency,  transient)  or  continuous  incident  wave  (single  fre¬ 
quency,  time  harmonic  steady  state) 

•  Numerical  grid  generation  —  structured  multizone  grid  or  unstructured  grid 

Complex  geometry  with  layered  radar  absorbing  media 

•  Lossy  or  lossless  material  properties 

Frequency  and  time  dependent  properties 
Thin  structures  (resistive  card,  lossy  pciint) 

•  Vector/parallel  code  mchitecture  —  2  GFLOPS  demonstrated  on  the  Cray-YMP 
with  8  processors,  and  10  GFLOPS  on  the  Cray-C90  with  16  processors.  Scalable 
performance  demonstrated  on  both  the  nCUBE  and  the  Intel  Paragon. 

Received  the  1990  CRAY  GigaBop  Performance  Award 

Received  the  1993  Computerworld  Smithsonian  Award 

•  Pre-  and  post-processor  graphics/animation 

•  Application  to  scattering  (RCS),  radiation  (antenna),  EMP/EMI/EMC,  and  bioelec¬ 
tromagnetics  problems 

•  Ideal  for  CFD/CEM  optimization  studies. 

The  CEM  code  has  been  extensively  tested  for  the  following  geometries. 

1)  Canonicail  objects  such  as  spheres,  cylinders,  ogives,  thin  rods,  cones,  airfoils,  and  a 
circular  disc 

2)  Almond  shaped  target 

3)  Inlets  of  various  shapes  (square,  circular,  curved,  •  •  •)  including  the  presence  of  infinite 
ground  plane 

4)  Flat  plates  of  various  planforms 

5)  Double  sphere 

6)  Complete  wing  geometries  with  layers 
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7)  Finned  projectile  and  cone-cylinder  combinations 

8)  Scattering  from  ship-like  targets 

9)  Complete  fighter  targets 

Some  sample  results  axe  shown  here  to  illustrate  the  present  capability.  Figure  6.1 
shows  an  unstructured  grid  set  up  for  an  almond.  Figure  6.2  shows  an  unstructured  grid- 
based  CEM  computation  for  a  complete  fighter.  Figure  6.3  shows  a  bioelectromagnetic 
computation  involving  the  aborption  of  microwave  radiation  by  a  human  used  in  the  hy¬ 
perthermia  treatment  of  cancer.  Figure  6.4  shows  resiilts  for  a  monopole  antenna.  Figure 
6.5  shows  results  for  a  two-  dimensional  photonic  band  gap  structure  used  as  a  filtering 
device.  Currently,  this  CEM  technology  is  being  applied  to  study  the  interference  of  pylon 
mountings  in  the  field  test  RCS  measurements  of  low  observable  platforms  of  interest  to 
RATSCAT  shown  in  Figure  6.6. 
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Figure  6.1  Unstructured  gridding  for  an  almond 


surface  gridding 


Figure  6.2  Unstructured  grid— based  CEM  for  a  complete  fighter 

35 


surface  electric 


Figure  6.3  Bioelectromagnetics  application  of  CEM 
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Photonic  Band 


Figure  6.5  Photonic  band  gap  periodic  structure 
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<  GHz: 


Figure  6.6  Application  of  CEM  for  studying  interference  of  pylon  mountings  in  RCS  measurements 
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7.0  MASSIVELY  PARALLEL  COMPUTING 


Parallel  Implementation 

With  the  emergence  of  massively  parallel  computing  architectures  with  potential  for 
terafiops  performance,  any  code  development  activity  must  effectively  utilize  the  computer 
architecture  in  achieving  the  proper  load  balance  with  minimiun  internodal  data  commu¬ 
nication. 

The  structured  finite-volume  code  was  originally  developed  and  optimized  for  vector 
computer  architectures.  The  implementation  of  the  code  on  a  distributed  memory  parallel 
architecture  was  accomplished  by  re-using  much  of  the  original  vector  code.  Additional 
coding  was  added  for  handling  inter-processor  communication  and  other  functions  unique 
to  the  parallel  implementation. 

Parallelization  Strategy 

For  the  structured  formulation  of  the  finite-volume  code  the  computational  domain 
surrounding  the  target  geometry  is  composed  of  3-dimensional  6  sided  voliimes  of  grid 
points  called  zones  or  blocks.  Each  side  or  face  of  a  zone  either  connects  to  another  zone  or 
has  a  boundary  condition  defined  on  that  face  (perfect  conducting  surface,  outer  boundary, 
etc.).  The  parallel  algorithm  takes  advantage  of  this  multi-zonal  gridding  capability  in 
order  to  divide  work  among  processors.  The  vzirious  zones  are  grouped  onto  processors, 
with  each  processor  obtaining  a  solution  for  the  cells  within  its  own  local  set  of  zones. 

Communication  Requirements 

The  solution  procedure  does  not^  allow  for  processors  to  proceed  completely  asyn¬ 
chronously.  Solving  for  cells  on  zone  faces  that  are  connected  to  other  zones  requires 
information  from  within  the  adjacent  zone.  This  information  may  be  available  locally  if 
the  adjacent  zone  resides  on  the  same  processor,  or  message  passing  may  be  required  if 
the  adjacent  zone  resides  on  another  processor.  This  boundary  update  message  passing  or 
flux  transfer  message  passing  is  done  twice  per  solution  time  step  and  forms  the  bulk  of 
the  parallel  code’s  message  peissing  requirements. 

Load  Balancing 

Load  balancing  is  achieved  by  mapping  zones  onto  processors.  Perfect  load  balancing 
requires  that  each  processor  have  the  same  number  of  zones,  each  containing  the  same 
number  of  grid  points  and  equal  numbers  and  types  of  boundary  condition  cells.  Simple 
geometries  may  usually  be  zoned  in  such  a  manner  as  to  obtain  perfect  load  balancing. 
For  complex  geometries  perfect  load  balancing  is  much  more  difficult,  but  adequate  load 
balancing  may  usually  be  obtained  by  mapping  a  close  to  equal  number  of  grid  points  onto 
each  processor. 

Scalability  Results 

Validation  and  timing  studies  have  been  performed  on  a  512-node  nCUBE  and  a 
208-node  Intel  Paragon.  Currently  the  code  shows  good  scalability  on  evenly  balanced 
test  cases.  These  cases  typically  had  simple  gridding  requirements  and  a  straight  forward 
domain  decomposition.  The  results  show  that  inter- processor  communication  due  to  flux 
transfer  never  becomes  a  dominant  time  factor  even  on  problems  with  large  numbers  of 
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grid  points  run  on  many  processors.  The  sphere  test  case  illustrated  in  Figure  7.1  shows 
how  problem  size  and  number  of  processors  can  be  increased  while  solution  time  remains 
level.  Perfectly  conducting  sphere  grids  were  nm  on  6,  24,  and  96  processors  of  the  Intel 
Paragon.  The  number  of  grid  points  per  processor  remained  constant  at  approximately 
60,000  resulting  in  total  grid  sizes  of  approximately  0.35, 1.4,  and  5.7  million  grid  points  for 
the  three  cases.  Since  increasing  the  number  of  processors  results  in  an  increased  number 
of  zonal  interfaces,  flux  message  passing  requirements  increase  throughout  the  system. 
Despite  this  increase  in  required  message  passing,  communication  times  did  not  change 
appreciably. 

Complex  problems  such  as  full  scale  fighter  geometries  also  show  encouraging  results. 
Figure  7.2  shows  timing  results  and  zoning  for  the  VFY218  fighter  gridded  for  a  frequency 
of  500MHz  with  a  10  point  per  wavelength  resolution.  A  total  of  58  zones  and  2.2  million 
grid  points  were  required.  The  grid  was  run  on  an  Intel  Paragon  using  28,  61,  and  128 
processors.  Preliminary  timing  data  reveals  that  communication  overhead  remains  at 
between  1  and  2.5  percent  of  the  total  solve  time  and  that  solution  speedup  occurs  as  the 
problem  is  distributed  over  more  processors.  Speedup  may  be  improved  by  addressing  load 
balancing  issues  arising  from  complex  zoning  arrangements. 
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Figure  7.1  Scaling  of  computer  resources  with  problem  size  for  a  sphere 


Parallel  CEM  for  a  Full  Fighter 


Figure  7.2  Scaling  of  computer  resources  with  problem  size  for  a  fighter 


8.  FURTHER  DEVELOPMENT  OF  THE  UNSTRUCTURED-GRID  CEM  CODE 


Based  on  the  preliminary  results  obtained  so  far  from  the  unstructured— grid  CEM 
code,  many  issues  need  further  study  to  matmre  this  technology  comparable  to  the  struc¬ 
tured  grid  RCS3D  code.  Some  of  them  axe, 

1.  The  process  of  constructing  a  polynomizd  basis  function  within  each  cell  for  the  electric 
and  magnetic  fields  requires  further  study  for  determining  the  accuracy.  Current 
procedure  involves  the  use  of  Riemann  fluxes  at  cell  interfaces  and  greatly  simplifies 
the  process  of  evaluating  the  slope  information.  For  arbitrary  cell  arrangements  this 
procedure  requires  further  work,  especially  at  boundary  faces. 

2.  We  need  to  look  at  higher  order  basis  functions  greater  than  second  order  accuracy. 
There  is  a  computational  penalty  that  goes  with  the  evaluation  of  higher  order  poly¬ 
nomials.  We  need  to  understand  the  balance  between  computational  efficiency  and 
accuracy. 

3.  Current  procedure  does  not  solve  the  Maxwell  equations  at  a  metal  boundary  (option 
10  type  procedure  in  RCS3D)  and  updates  the  electric  and  magnetic  fields  only  at 
cell  centroids  as  a  cell  averaged  value.  We  need  to  look  at  ways  to  improve  the 
accuracy,  especially  for  traveling  waves,  by  solving  Maxwell’s  equations  right  at  metal 
boundaries  within  the  unstructured  grid  approach  (option  1  in  RCS3D). 

4.  The  current  code  can  have  only  one  type  of  a  cell  arrangemnt  in  the  entire  computa¬ 
tional  domain  (hexahedron,  prism,  or  tetrahedron).  For  geometric  flexibility  we  need 
to  upgrade  the  pilot  code  to  inclilde  combinations  of  different  cell  shapes. 

5.  The  present  algorithm  reduces  to  Lax— Wendroff  for  hexahedral  cell  shapes.  This 
allows  one  to  construct  a  hybrid  structured/unstructured  grid  solver  to  exploit  the 
virtues  of  both  approaches. 

6.  RCS3D  has  many  types  of  boundary  conditions  such  as  impedance  layers,  resistive 
cards,  and  can  handle  frequency  (dispersive)  and  time  dependent  (active)  material 
properties.  The  unstructured  CEM  needs  further  development  to  include  all  these 
features. 

7.  Pre  and  post  processing  for  unstructured  arrangement  needs  further  work.  The  FAST 
capability  developed  by  NASA  Ames  Research  Center  is  rather  limited  and  cannot 
handle  hexahedron,  prism,  and  tetrahedron  combinations.  Also,  it  cannot  handle 
plotting  of  electric  and  magnetic  fields  at  cell  centroids  while  the  grid  is  specified  at 
nodes. 

8.  The  use  of  tetrahedra  can  result  in  small  cell  volumes  (compared  to  a  similar  hexa¬ 
hedron).  This  can  cause  the  time  step  size  to  be  small  (maybe  a  factor  of  four  less 
compared  to  hexahedra  cells).  We  need  to  look  at  pointwise  implicit  schemes  to  allow 
larger  time  steps  while  maintaining  stability. 

9.  The  current  “advancing  front”  grid  generator  needs  further  work  to  include  constraints 
on  minimal  allowable  cell  volumes.  In  the  present  approach,  arbitrarily  small  volumes 
can  result  at  a  distance  from  the  scattering  object  causing  small  time  steps. 
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10.  We  need  further  study  to  understand  the  grid  resolution  requirements  in  the  unstruc¬ 
tured  tetrahedral  arrangements  (how  many  cells  per  wavelength?). 

11.  Unstructured  algorithms  are  inherently  not  well  suited  for  optimum  vectorization 
because  of  indirect  addressing.  For  Cray-hke  machines  (coarse  grain  vector /parallel), 
we  need  to  understand  the  proper  data  and  code  structure  to  acieve  the  type  of 
MFLOP  ratings  demonstrated  by  RCS3D.  The  current  version  runs  about  three  times 
slower  than  RCS3D. 

12.  The  real  power  of  an  unstructured  grid  approach  is  in  using  MIMD  architectures. 
We  need  to  understand  how  to  achieve  proper  load  balancing  and  choices  for  do¬ 
main  decomposition.  This  is  a  crucial  research  topic  for  demonstrating  large  RCS 
computations  such  as  a  complete  fighter  at  giga  Hertz  frequencies. 
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Al.  Treatment  for  Impedance  Layer 


I.  INTRODUCTION 

The  electromagnetic  signatures  of  perfectly  conducting  objects  may  be  reduced  by 
coating  their  surface  with  a  thin  layer  of  lossy  dielectric,  absorbing  material  such  as  paint. 
This,  together  with  the  shape  factor,  axe  the  basic  tools  for  developing  low  observable 
structures.  In  computing  the  radax  cross  section  (RCS)  of  these  objects,  if  the  brute 
force  approach  is  attempted,  the  problem  is  often  too  large  to  be  handled,  since  effective 
coatings  usually  have  high  complex  permeabilities  and  permittivities.  Therefore,  some 
form  of  approximate  boundary  conditions^®  have  to  be  used  to  replace  the  dielectric  layer 
by  an  impedance  surface  at  the  interface  between  the  dielectric  and  free  space. 

The  impedance  boundary  condition,  however,  is  a  frequency  domain  concept.  In 
this  appendix,  it  is  shown  how  this  concept  can  be  extended  to  the  time  domain  and 
implemented  for  the  finite-volume  time— domain  (FVTD)  method. 

II.  FORMULATION 

The  most  basic  problem  to  which  the  impedance  boundary  condition  can  be  applied 
is  that  of  a  plane  wave  incident  upon  a  material  half— space.  Provided  that  the  magnitude 
of  the  complex  refractive  index  of  material  is  large,  i.e.,  \N\  »  1  where  N  = 
with  fj.,  e,  /xo,  and  €o  being  the  complex  permeabilities  and  permittivities  of  the  material 
half-space  and  free  space,  respectively,  then  the  material  interface  may  be  replaced  by  an 
impedance  surface  located  at  the  interface  on  which  the  electric  and  magnetic  fields  satisfy 
the  following  condition,  known  as  the  impedance  boundary  condition: 

n  X  (n  X  jF)  =  -r]h  x  H  (^Ul) 

where  rj  is  the  characteristic  impedance  of  the  material  medium  and  n  is  a  unit  vector 
normal  to  the  interface  in  the  free  space.  This  impedance  boundary  condition  remains  valid 
when  the  interface  is  a  curved  surface  provided  the  additional  constraint  |Im(iV)|/i:oP  >>  1 
is  imposed,  p  is  the  minimum  radius  of  the  curvature. 

This  basic  concept  may  be  simply  extended  to  a  multilayered  medium.  Only  the 
simple  case  of  a  lossy  material  coating  over  a  conducting  body  will  be  discussed  here. 
From  a  simple  transmission  line  anzdogy,  the  coated  object  may  be  replaced  by  a  surface 
whose  impedance  is  given  by: 

T]  = —iZ  taji{kody/fIe)  (A1.2) 

where  Z  is  the  characteristic  impedance  of  the  material  coating,  d  is  its  thickness,  ko  is  the 
propagation  constant  in  free  space,  and  a  time  dependence  exp(—iut)  has  been  assumed. 

In  the  framework  of  the  FVTD  method,  Maxwell  s  equations  in  the  partial  differential 
form  in  a  body-fitted  coordinate  system  are  numerically  integrated  in  the  time  and  over 
the  faces  of  each  grid  cell  to  compute  the  field  values  inside  that  cell,  as  explained  in  the 
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previous  sections.  When  one  or  more  of  the  faces  of  a  cell  are  impedance  surfaces,  then 
a  special  treatment  is  necessary.  The  tangential  components  of  the  electric  and  magnetic 
fields  of  each  cell  used  in  the  time  integration  are  derived  from  jump  conditions  usually 
defined  in  the  time  domain.  In  the  present  case,  however,  it  advantageous  to  start  with 
jump  conditions  in  the  frequency  domain  and  derive  the  expressions  for  the  tangential  fields 
on  the  impedance  faces  of  the  cells  involved.  Then,  the  time-domain  tangential  fields  are 
derived  from  their  frequency-domain  counterparts  through  Fourier  transformations. 

Figure  Al.l  shows  a  cell  in  the  computational  grid.  The  tangential  fields  for  this  cell 
with  one  of  its  ^-constant  faces  being  an  impedance  surface  will  be  derived.  Derivation  for 
rj-  and  ^-constant  faces  are  identical.  The  jump  condition  across  the  positive  characteristic 
as  depicted  in  Figure  Al.l  may  be  expressed  as: 

(^H+  -  H*^  =  C  X  (e+  -  (A1.3a) 

where  the  field  quantities  bearing  superscript  ‘+’  are  defined  at  the  centroid,  those  with 
asterisk  on  the  faces  of  the  cell,  and  is  the  characteristic  impedance  of  the  material 
filling  the  cell.  Equation  (A1.3a)  basically  relates  the  discontinuities  of  the  electric  and 
magnetic  fields  across  the  eigenvalues  A"*"  =  c|^|  of  Maxwell’s  equations  where  c  is  the 
speed  of  light  in  the  medium.  For  both  the  scattered  and  total  fields,  (Al.3a)  assumes  the 
same  form.  However,  it  is  meant  to  present  scattered  fields  here.  In  terms  of  the  scattered 
fields,  (Al.l)  may  be  written  as  follows: 


i  X  {E*  +  iol  =  X  {H*  +  H‘) 


(.41.36) 


where  rj  is  given  by  (A1.2).  From  (A1.3)  (  x  H*  is  obtained. 


-  n  +  z* 


(yll.4) 


where  ^  x  It  is  advantageous,  however,  to  find  |  x  E*  from  a  jump 

condition  similar  to  (A1.3a)  in  the  time  domain.  This  will  be  done  later.  Inverse  Fourier 
transformation  of  (A1.4)  results  in  the  desired  expression  for  the  tangential  magnetic  field 
on  the  impedance  face  of  the  cell  in  the  time  domain  which  is  as  follows: 


I  X  H\t)  =  [  F{i-  T)i  X  VV+(r)dr 

Jo 

+  f  F(f  -  r)|  X  >V‘(r)dr  -  I  X  ?i‘(t) 

Jo 


(Al.o) 
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where: 


m 

w 

Re[Tj] 

Im[rj] 

D(u) 

Ni(u) 

N2{u) 

a 


'{z^)  ’ 


W)  ’ 

jNi{u)  +  6N2(u) 

Diuj) 

1  +  tan^(kodQ)  t&nh^{kod^)  , 
taii(koda)/  cosh^(fcorf/?)  , 
iaiih(kodl3)/  cos^{koda)  , 

=  ±yj  -(i2e[/xe]  +  l^e|)  , 

P  =  ±y^(--Re[/ie]  +  |/xe|)  , 


7 

<5 


± 

± 


Fig.  Al.l.  A  Cell  with  Impedance  Face. 

The  plus  signs  apply  to  a  wave  traveling  in  the  positive  coordinate  direction,  while 
the  minus  signs  to  a  wave  traveling  in  the  negative  coordinate  direction. 

For  a  time-harmonic  or  narrow-band  plane  wave  excitation,  provided  the  steady  state 
solution  is  desired,  a  fictitious  dynamics  for  rj{uj)  will  be  constructed  whose  inverse  Fourier 
transform  can  be  evaluated  in  a  closed  form.  The  following  two  cases  will  be  considered 
here: 
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1)  ^^<0 


In  this  case  we  let: 


rj{uj)  =  r  —  iux 


(^1.6) 


where  r  and  x  axe  constant  over  the  frequency  band  of  interest.  The  straightforward 
application  of  an  inverse  Fourier  transform  results  in  the  desired  function  in  the  time 
domain 

F{t)  =  ae"^‘  (A1.7a) 


where: 


1 

a  =  — 

X 

r  +  Z+ 


2) 

r){uj)  in  this  case  will  be  considered  to  be: 


(Al.lb) 


(Al.lc) 


ri{uj)  =  r'  —  i —  . 

U) 

The  inverse  Fourier  transform  of  (A1.8)  results  in  the  following  function: 


where; 


F{t)  —  cS{t)  +  ae' 


Z+  +  r' 


(Al.S) 


(A1.9a) 


(A1.96) 


+  r' 


r'  +  Z+ 


(A1.9c) 


(AIM) 


For  notation  simplification  we  use  (Al.9a)  for  both  cases  with  c  =  0  for  Case  1. 
A  plane,  time-harmonic  incident  field  may  be  represented  by: 

£i{t)  =  9  cos{k  ■  R  —  u't) 


(Al.lOa) 


'Hi(t)  =  cos{k  ■  R  —  u:t)  . 

£j  ' 


From  (Al.lO),  W'{t)  may  be  constructed: 


W\t)  =  Z+H\t)  -ix  E'{t) 
=  Cj  cos{K  ■  R  —  ijjt) 


(A1.106) 


(Xl.ll) 
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where  =  ^  + 1  x  (9. 

Equations  (A1.9),  (Al.lO),  and  (Al.ll)  axe  substituted  in  (A1.5)  to  obtain  the  mag¬ 
netic  field  fiux  term  in  the  time  domain  for  the  impedance  surface: 


( X  H*{t)  =  4  X  +  4  X  W\t) 

+  a  f  X  W*(r)dr  (_41  12) 

Jo 

+  a  f  X  W\T)dr  -  |  x  W‘(t) 

Jo 

The  first  integral  in  the  right  hand  side  of  (A1.12)  will  result  in  a  recursive  algorithm.  The 
details  follow.  Let:  ^ 

f{t)  =  /  X  W+{T)dT  .  (A1.13) 

Jo 

f{t)  may  be  regarded  as  a  formal  solution  to  the  following  first  order  differential  equation: 

^  +  bf  =  ixW+{T)  .  (A1.14) 

at 

It  is  found  that  the  following  implicit  algorithm  can  provide  accurate  solution  to  T  for  a 
large  range  of  values  of  b:  ^ 

f{t  +  At)  -  f{t)  b  j ^  [I  X  W-^{t  -h  At)  +  I  X  #+(t)j  .  (A1.15) 

At  2  ^  ^ 

This  will  result  in  the  following  recursive  formula  for  computing  T (t)  for  the  discrete  time 
variable: 


2-bAt^ 

2-f6At 


At 

2  +  bAt 


{i  X  w+y  +  {i  X  w+y  ^ 


(A1.16) 


where:  ^ 

T„  =  T{nAt) 

To  =  0 


(I  X  W^y  =  lx  W+{nAt) 

The  second  integral  in  (A1.12)  can  be  evaluated  in  a  closed  form.  Let: 

U{t)  =  f  X  W’{T)dT 

Jo 

=  ixW  [  cos(jb  •  R  -  WT)dT 

Jo 


(.41.17) 
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Then  it  jnay  be  shown  that: 


u{t)  =  ixw- 


cos(fc  •  R  —  u>t) 


\r=—<j> 


(>11.18) 


I  J 


Expressions  (Al.16)  and  (Al.18)  are  substituted  in  (A1.12)  to  obtain  the  desired  magnetic 
field  fiux: 

ixn*{t)  =  cixW^{t)-\-cixWi(t)  +  af{t)  +  aU(t)-ixH*(t)  .  (A1.19) 


The  results  in  (A1.19)  may  be  presented  using  the  discretized  time  variable  n: 


(I  X  n*y  =  c(ix  W+y  +  C  (e  X  W'y  +  af„  +  aUn  -  (l  X  h)’  (Al.20) 

where  tjn  =  U(nAt). 

The  electric  field  flux,  ^  x  £*,  may  be  obtained  directly  from  the  cross  product  of  the 
unit  vector  ^  and  the  following  jump  condition  in  the  time  domain: 

r  -  £+  =  I  X  [z+(W*  -  '«+)]  .  (A1.21) 
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A2.  Treatment  for  Anisotropic  Media 


Radar  absorbing  materials,  such  as  filled  honeycomb  structures,  may  exhibit  consid¬ 
erable  anisotropy  in  their  response  to  electromagnetic  fields.  Treating  these  structures  as 
continuous  media  requires  tensor  definitions  for  e  and  n-  Thus,  the  interface  fluxes  derived 
earlier  for  isotropic  media  need  to  be  generalized.  We  start  from  the  conservation-law 
form  of  Maxwell’s  equations  in  curvilinear  coordinates  equation  (32): 


dt 


ixE/J 

-ixHU)  dr)\ 


ffxElJ  \,±( 
-ifxHij)  dC\ 


C^E/J  \ 
-CxH/j) 


where  Q  =  (B,D)  and  S  =  (0,-7).  As  in  the  isotropic  case,  we  shall  solve  the  one- 
dimensional  Riemann  problem  associated  with  one  coordinate  direction,  and  we  consider 
that  the  finite— volume  cell  is  bounded  by  surfaces  of  constant  and 

Writing  E  =  e-^D  and  H  =  one  finds,  for  example,  that  the  ^-derivative  term 

above  takes  the  form: 


(xE/J 

ixHU)  di 


X  e-' 

[-(xfi  ^0 


9.\  ^  —  [ A— 
j  ~  d(  rj 


(A2.1) 


where  the  position  dependence  of  e  and  fj,  is  accounted  for  in  the  6x6  matrix  A.  In  our 
finite-volume  formalism,  e  and  fi  areT,-aken  to  be  constant  inside  each  cell,  and  the  fields 
B  and  D  are  evaluated  at  the  cell  centroid. 


Waves  that  propagate  along  the  f  direction  correspond  to  eigenvectors  of  A.  As  we 
shall  see,  these  eigenvectors  are  pairs  (^BaiEo^  having  both  Ba  and  Ea  lying  in  the 

plane  perpendicular  to  f  at  the  interface  between  neighboring  cells.  Because  the  material 
properties  may  vary  from  cell  to  cell,  these  eigenvectors  in  general  also  change  across  the 
interface. 

In  terras  of  the  3  x  3  matrix  5  =  x  x  ,  the  right  eigenvectors  and 

eigenvalues  on  one  side  of  a  ^—interface  can  be  derived  in  the  form. 


Cl  =  -  (5ii  +  S22  +  *533)  ;  C2  =  -S12S21  -  Si^Szi  -  S23SZ2  +  S11S22  +  S22S33  +  5ii533  , 

2  _ _  (A2.2c) 

Zi  =  ^yiv2  '  e-H2)  +  {vi  •  ;  Z2  =  y/{vi  ■  I  {v2  ■  fi~'^V2)  •  (A2.2d) 
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The  physical  tinit  vectors  v\  and  are  perpendicular  to  and  they  are  determined  by 
the  requirements: 

Ul  •  p~^V2  =  V2  •  ^2  *  e“^Ul  =  Vl  •  €~^V2  =  0  .  (j42.3) 

The  eigenvectors  rs  and  re  do  not  correspond  to  propagating  waves  because  As  =  Ae  =  0, 
and  they  are  not  required  in  constructing  the  interface  flux. 

As  in  the  isotropic  case,  the  basic  relations  constraining  the  interface  flux  are  that 
the  difference  between  the  centroidal  values  of  B  and  D  and  their  values  at  the  interface 
must  lie  along  the  eigenvectors  (a  =  1  to  4).  Coupled  with  the  condition  that  the 
components  of  E  and  H  tangential  to  the  interface  must  be  the  same  on  either  side,  these 
relations  determine  the  following  expressions  for  the  interface  fluxes: 


-r}(a+  +  a-)  ■[&+  -R-  +  N  x^-  C(b+  +  b-)[H+  -  g"  +  iV  x  |] 
I  •  (a+  +  a“ )  X  (b+  +  b~) 


(A2.4a) 


ixH*  = 

Tj{Z+Z}a+i-Z^Z^a-)-[E+-E--M  x  i\  +  ({Z+Z+b++Z^Z^b-)-[E+  -  E'-M  x  |] 


i  •  (Z^Z^a+  +  X^Z^a-)  x  (Z+Z+h  +  ZiZ^h) 
where  the  vectors  a  and  b  are  given  by: 


(A2.46) 


^  _  Zi{vi  •  7?)vi  +  Z2{v2  ■  1)^2  .  ^  _  Zi{vi  •  C)r)l  +  Z2{v2 ' 0^2 

Z\Z2^  •  (vi  X  V2)  Z\Z2^  •  (vi  X  t'2) 


{A2.d) 


Sind  the  vectors  M  and  iV  axe: 

M  =  Z+Z+  [a^{H+  ■  77)  +  b+{H+  ■  0]  +  zrz^  +  b-{H  ■  O]  (A2.6a) 

N  =  [a+(£+  •  77)  +  b+{E+  •  0]  +  [a~{E-  ■  77)  +  r  (i"  •  o]  •  (^2.7a) 
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Appendix  A3.  Geometry  and  Gridding  Techniques 


This  appendix  contains  Sections  12-14  of  a  Navy  contract  report,  NAWCWPNS  TP 
8220,  prepared  by  the  Computational  Fluid  Dynamics  group  at  Rockwell  Science  Center. 
The  following  page  numbers  have  been  taken  directly  from  that  report. 
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NAWCWPNS  TP  8220 


12.0  GEOMETRY  FORMULATION  DETAILS 

The  geometry  formulation  in  UNIVERSE-series  unstructured  grid  codes  is  very 
flexible  and,  in  general,  can  handle  any  type  of  element  (conservation  cell).  In  partic- 
tdar,  we  have  included  three  types  of  cells  at  present.  These  are  1)  the  hexahedral  cell, 
2)  the  triangular  prism  ceU,  and  3)  the  tetrahedral  cell.  The  details  of  the  geometry 
treatment  for  these  three  cell  types  are  given  in  this  section. 


12-1  Local  Node  Numbers  and  Coordinates 

A  cell  type  is  defined  by  the  number  and  placement  of  its  vertices.  A  complete 
definition  of  a  cell  type  would  also  include,  in  addition  to  the  vertices,  any  other 
nodes  needed  to  define  the  variation  of  the  geometric  variables  i,y,r  in  the  cell  (see 
Fig.  4.2  on  higher-order  elements).  Taken  together,  the  vertex  nodes  and  other  nodes 
<-aTi  simply  be  referred  to  as  the  nodes  of  the  cell. 

We  assume  a  local  coordinate  system  for  each  cell  (^,t?,<7,  Fig.  4.1).  The  vertex 
node  numbers  are  defined  in  terms  of  an  ordered  set  of  integers  and  the  values  of  the 
local  coordinates  at  each  vertex  are  defined  for  each  cell  type  in  Figs.  12.1,  12.2  and 

12.3. 

12.2  Shape  Function 

The  geometric  “shape  function”  Ni{i,r),a)  was  defined  in  Eqs.  4.9-11  for  the 
three  cell  types.  Each  cell  node  is  associated  with  its  own  shape  function  which  is  a 
multidimensional  polynomial  in 

k=^0 

The  polynomial  coeffidents  can  easily  be  obtained  by  solving  the  K  equations  result- 
ing  from  setting 
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(12.2a) 

(12.2a) 


Vertex  nodes  and  coordinates  for  hexahedron 
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Figure  12.3  Vertex  nodes  and  coordinates  for  tetrahedron 
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meu  the  th«e.dime»sionel  geometry  polynomial  is  bnilt  by  forming  the  product  o  > 

t^oormorelower-dimensional  polynomials  (Eq.  4.10, Eq.  4.11), the  three.d.mens.onel 

shape  functions  can  also  be  constructed  from  the  appropriate  lower-d.mens.onal  shape 
functions. 


12.2.1  T.ncal  Faces  apd  Vertices 

Each  face  of  a  ceU  is  defined  by  which  vertices  belong  to  it  (in  terms  of  local 
node  numbers).  The  faces  and  their  corresponding  vertices  are  identified  for  each  of 
the  three  cell  types  in  Figs.  12.4, 12.5  and  12.6. 


12.2.2  C«»11-Face  Metrics  and  Normals 


The  method  used  to  compute  cell-face  normals  was  outlined  in  Section  4  in  the 
paragraphs  following  Eq.  4.11.  A  more  detailed  description  is  given  below. 

Given  the  geometry  polynomials  (Eq.  4.9,  Eq.  4.10,  or  Eq.  4.11),  the  tangents  to 
the  local  coordinate  directions  can  be  defined  by 


^  +  y^k  +  zi'l 

rf  =  X,jj  +  t/jjfc  + 

5  =  Xcj  +  y^k  Zffl 


(12.3) 


Neat  we  identify  two  vectors  in  the  plane  of  the  ceU  face  by  choosing  suitable  linear 
combinations  of  the  above. 


®^(U\  h2  <13^ 
V<21  <22  <23/ 


(12.4) 


The  transformation  matrix  in  the  above  can  be  defined  as  T 


/<n  <12  <13  \  (12.5) 

\t21  <22  <23/ 


and  is  defined  for  each  face  of  each  cell  type  in  Figs.  12.7,  12.8  and  12.9.  The  T 
matrix  elements  are  selected  appropriately  to  account  for  the  fact  that  range. 
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in  each  cell,  from  -1  to  +1  in  many  cases.  The  cross  product  of  Vi  and  V2  defines 
the  cell-face  normail. 

12.3  Cell-Face  Quadrature  Points 

At  every  cell  face,  the  integral  of  any  quantity  Q  over  that  face  can  be  replaced, 
up  to  the  desired  order  of  accmacy,  by  a  suitable  quadrature  formula 

f  I  Q(s,t)ds  dt  =  E  CiQisi,ti)  (12.6) 

•'»  quads. 

where  s,t  are  the  representative  running  orthogonal  coordinates  tangential  to  the 
face,  and  subscript  i  denotes  a  quadrature  point.  The  coefficient  C,-  is  the  “weight” 
assigned  to  the  i-th  quadrature  point.  For  the  three  cell  types  being  considered,  one 
encounters  only  two  types  of  cell  faces  —  “square”  or  “triangular  (in  terms  of  the 
running  coordinates  s  and  f,  as  well  as  the  local  coordinates  ^,T],cr).  Tables  12.1  and 
12.2  represent  four-point  Gaussian  quadrature  points  and  weights  for  such  “square 
and  “triangular”  cells.  This  leads  to  fourth-order  accuracy  in  the  evaluation  of  cell- 
face  integrals.  Figure  12.10  presents  these  quadrature  points  diagrammatically  (not 
necessarily  to  true  scale). 


s 

t 

C 

-0.577350269189626 

-0.577350269189626 

0.5 

-1-0.577350269189626 

-0.577350269189626 

0.5 

-0.577350269189626 

-1-0.577350269189626 

0.5 

-1-0.577350269189626 

+0.577350269189626 

0.5 

Table  12.1  Gaussian  quadrature  points  for  “square”  face 
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s 

t 

C 

-0.487831473351 

-0.689897894859 

0.318041443825 

-1-0.487831473351 

-0.689897894859 

0.318041443825 

-0.204988807440 

+0.289898127317 

0.181958556175 

-1-0.204988807440 

+0.289898127317 

0.181958556175 

Table  12.2  Gaussian  quadrature  points  for  “triangular  face 


12.3.1  Special  Cases 

While,  in  general,  the  four-p&int  Gaussian  quadrature  is  employed  for  three- 
dimensional  problems,  simpler  formulae  (with  fewer  quadrature  points)  can  be  used 
sometimes  without  any  loss  of  accuracy.  For  example,  with  hexahedral  cells,  the 
midpoint  rule  is  sufficient  for  one-dimensional  problems  in  the  one  significant  direction 
being  considered.  For  two-dimensional  problems,  a  two-point  quadrature  is  sufficient 
for  faces  spanning  the  third  direction  (the  two  points  are  along  the  line  that  bifurcates 

the  cell  in  the  third  direction). 


12.4  Cell  edges 


The  need  to  determine  the  intersection  between  interior  surface  geometry  ele¬ 
ments  and  elements  that  define  a  cutting  plane  or  surface  for  postprocessing  purposes 
required  a  numbering  system  for  cell  edges  to  be  added  to  the  book-keeping  frame¬ 
work.  Figures  12.11,  12.12  and  12.13  illustrate  the  vertex  nodes  and  cell  edges  for 
the  three  types  of  cells  considered  in  this  report. 
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Calnilation  of  Cell  Volume 


The  cell  volume  is  related  to  Eqs.  4.18  and  4.19 

-///.  V  ■{xj  +  yk  +  zl)dV 


(12.7) 


corresponding  to  t  =  0  in  Eq.  4.18.  The  volume  integration  can  be  replaced  by  surface 

integration.  r  t  r 

y  =  /  /  /  iV‘Xo)dV 

J  J  Jv  (12.8) 

=  J  J{XQ-n)dS 

The  surface  integral  can,  in  turn,  be  converted  to  a  quadrature  formula  based  on 
Eq.  12.6. 

12-6  Strong-Conservation-Law  Form  in  General  Coordinates 


In  Section  5,  the  strong-conservation-law  form  was  introduced  for  hyperbolic 
systems  of  conservation  laws  in  Cartesian  coordinates.  In  future  editions  of  this 
report,  the  invariance  of  that  form  under  coordinate  transformations  will  be  discussed. 

12.7  Integral  and  Strong-Conservation-Law  Forms 


There  is  a  clear  and  useful  geometric  analogy  between  the  integral  form  and 
the  strong-conservation-law  form  of  the  equations.  This  will  be  elucidated  in  future 
editions  of  this  report. 

12.8  Preservation  of  Uniform  Flow 


Physically,  imiform  flow  remains  uniform.  Numerically,  this  simple  fact  may 
not  be  true.  In  other  words,  when  a  numerical  method  is  applied  to  update  the 
solution  from  time  step  to  time  step  with  uniform  flow  as  the  starting  solution,  on 
nonuniform  grids  only  a  careful  construction  of  the  algorithm  will  guarantee  that  the 
numerical  solution  does  not  introduce  spurious  disturbances  into  uniform  flow.  All 
that  is  required  is  that  formulae  of  high  enough  order  of  accuracy  are  employed  in  the 
boundary  quadrature  procedure.  These  issues  will  be  considered  in  detail  in  future 
editions  of  this  report. 
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Figure  12.6  Faces  and  their  vertices  for  tetrahedron 


NAWCWPNSTP8220 


2 


Figure  12.8 


Faces  and  their  transformation  matrix  for  tnangular  prism  continued 
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Figure  12.8  Faces  and  their  transformation  matrix  for  triangular  prism  —  continued 
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Face  3 


Figure  12.9  Faces  and  their  transformation  matrix  for  tetrahedron  —  continued 
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Quadrature 
points  for 
4-vertex  face 


Quadrature 
points  for 
3-vertex  face 


Figure  12.10 


Gaussian  quadrature  points  for  “square”  and  “triangular”  faces 
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Figtire  12.11  Vertex  nodes  and  cell  edges  for  hexahedron 
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I 


Figure  12.12  Vertex  nodes  and  cell  edges  for  triangular  pnsm 
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13.0  MESH  GENERATION 

Th.  computer  code  known  as  TRIM3D  in  Phase  1 

L  unstructured  surface  srid,  viscous  triansular  P-  e 

grid  and  setup  of  boundary  conditions  are  generated  mteracUvely  rn  UNISG.  ^ 
Sd  generator  for  tetrahedral  volume  ceUs  is  UNIVG.  The  flow  solver  rs  UNIV.  The 
postprocessing  mode  of  UNIV  is  referred  to  as  UNIVP. 

The  computer  codes  UNISG  and  UNIVG,  based  on  the  two  layer  rnethod  to 
generate  3-dimensional  unstructured  grids  for  viscous  and  inviscd  flows  ha,«  been 
Leloped  at  RockweU  Science  Center.  The  boundary  of  the  computaf  onal  domarn 
is  divided  into  several  patches.  The  geometry  of  eadr  patch  is  d«cnbed  ’’X 
a  sufficient  number  of  non-intersecting  lines  on  the  pate  .  M  me,  1  , 

described  in  terms  of  a  sufficient  number  of  points  on  the  hne  (Fig.  13.1).  The  p 
^fLation  can  also  be  provided  in  the  form  of  IGES  NURBS  surfaces.  Flom  t^i 
information  a  parametric  representation  (s,  t)  of  the  surfa«  is  developed  where  imd 
i  are  the  local  running  coordinates  in  the  plane  of  the  surface.  For  convenience,^ 
parameters  s  and  t  are  chosen  to  take  values  between  0  and  1,  where  s  =  0  s  -  1. 

t  =  0andl  =  1  form  the  segmented  boundaryofthe  patch.  No  esare  en  is 

on  these  four  boundary  lines  satisfying  the  user-specified  clustering  -A— 
the  code,  the  first  and  last  hues  of  input  patch  information  corr^pon  o 
t=l  respectively.  Presently,  the  code  requires  that  the  number  of  nodes  at  le« 
«iy  one  pair  of  opposite  sides  (s  =  0  and  s  =  1  or  1  =  0  and  i  =  1)  be  equal.  T 
restriction  will  be  removed  in  the  future. 

Internal  lines  are  generated  by  connecting  corresponding  points  on  a  pmr  of 
opponte  sides.  Each  internal  line  is  then  divided  into  line  segments.  “  '^men 
length  at  the  ends  of  an  internal  line  is  determined  by  the  average  of  he  length  of 
irete  segments  on  the  boundary  of  the  patch  that  intersec.  the  ^ven  mt.n^ 
line.  The  length  of  the  remaining  line  segments  on  an  internal  hne  ,s  obtained  y 
smoothly  blending  the  lengths  of  line  segments  at  the  ends.  Triangular  »s^am 
then  generated  by  connecting  the  line  segments  accor  mg  to  e  recipe 

Fig.  13.2. 
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THe  next  step  is  to  generate  the  boxmdary  that  separates  the  viscous  and  the 
inviscid  regions.  For  this  purpose  points  are  generated  along  each  patch  boundary 
at  a  specified  distance  along  the  normal  to  the  patch.  These  points  are  connected 
to  form  the  boundary  of  the  viscous  region.  This  boundary  is  then  smoothed  using 
an  elliptic  smoothing  procedure  and  transfinite  interpolation  is  used  to  generate  a 
smooth  surface  that  is  enclosed  by  it.  Projection  of  triangles  on  the  patch  surface 
onto  the  newly  generated  smooth  surface  is  performed.  Corresponding  triangles  on 
the  two  surfaces  are  connected  to  form  triangular  prisms.  These  prisms  are  first 
broken  into  smaller  prisms  and  then  divided  to  form  tetrahedra.  This  process  is 
carried  out  in  such  a  way  that  the  grid  clustering  required  to  resolve  small  scales 
near  the  body  surface  is  achieved.  All  the  connections  between  patches,  direction 
of  volume  cell  generation,  and  grid  clustering  information  are  interactively  setup  in 
UNISG  and  automatically  output  in  a  form  compatible  directly  with  UNIVG. 

In  UNIVG,  the  “advancing  front”  method  (Ref.  45)  is  used  to  fill  up  the  com¬ 
putational  domain  with  tetrahedra  starting  from  triangular  elements  on  the  viscous 
boundary.  This  technique  requires  specification  of  an  initial  surface  divided  into  tri¬ 
angular  elements.  The  boundary  of  the  computational  domain  is  used  as  the  intial 
surface  of  the  advancing  front.  The  smallest  element  on  the  surface  is  replaced  by 
the  three  sides  of  a  new  tetrediedron  constructed  with  that  element  as  the  base.  The 
fo\irth  node  of  the  tetrahedron  is  obtained  either  from  one  of  the  existing  nodes  or 
by  creating  a  new  one  (Fig.  13.3). 

Grid  generation  in  the  inviscid  region  consists  of  the  following  steps; 

(1)  Define  the  boundatry  patches  (including  body  surface  and  outer  boundary). 

(2)  Triangular  elements  are  then  generated  on  each  patch  based  on  user  supplied 
information  regarding  the  number  of  points  and  mesh  streching. 

(3)  These  triangular  elements  are  used  ais  the  initial  surface  of  the  ad^'ancing  front. 

(4)  Tetrahedral  cells  are  then  constructed  with  the  triangular  elements  as  bases.  In 
order  to  avoid  large  cells  intersecting  with  small  cells,  the  smallest  triangular 
element  from  the  list  of  faces  is  always  used. 
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(5)SdectororeateapcinttofonnaUtrahedralcdlfron.thebasetria.gu^^ 

(ABC,  in  Fig.  13.4). 

a )  The  point  D  can  belong  to  a  neighboring  triangulax  element  (ACD)  ^  W 
■  ae  it  can  generate  a  tetrahedron  satisfying  certain  built-rn  constrarnts,  (T 
angles  between  the  base  triangular  element  ABC  and  newly  o™'  ^ 
gles  ABD,  BCD,  CAD  should  be  less  than  135  degrees  and  larger  than 

degrees.)  See  Fig.  13.4a  (top  figure). 

b  )  The  point  D  can  he  selected  from  an  exbting  grid  or  from  a  hst  of  pomts 
suppUed  for  the  purpose.  See  Fig.  13.4b  (middle  figure). 

c  )  Otherwise,  a  new  point  D  located  on  the  line  which  is  normal  to  the  base 
triangular  element  at  its  center  can  be  created  to  form  the  new  tetrahedron. 

See  Fig.  13.4c  (lower  figure). 

(6)  Check  whether  the  newly  formed  cell  intersects  with  any  already  existing  cells. 

(7)  If  no  suitable  point  can  be  found,  some  old  cells  have  to  be  removed,and  repeat 
step  (5). 

(8)  The  new  faces  ABD  and  BCD  from  (5-a)  or  ABD,  BCD,  CAD  «b)  »d 
(5.C),  the  new  cell  ABCD  and  new  point  D  from  (S-c)  are  all  added 

respective  lists  (face,  cell  and  point). 

(9)  Delete  the  old  face  (ABC)  from  the  list. 

(10)  Check  if  there  are  any  triangular  elements  remmning  in  the  front  .  If  yes,  go 
back  to  step  4. 

.  r.  1  /.4c.  nf  TTNTVG  are  1^  data  structure  and  2)  the 

The  two  most  important  aspects  of  U  Ml  VO  ) 

lope  to  check  whether  any  two  tetrahedra  intersect  with  eadr  other.  These  aspe 
are  disscussed  in  the  following  section. 

nafH  Structure 


In  order  to  find  the 


smallest  triangular  element,  points  in  the  neighborhood  and 
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neighboring  triangular  elements  efficiently,  the  data  structure  for  the  faces,  and  the 
grid  points  is  arranged  as  follows: 

13.1.1  List  for  The  Front  Face 

A  binary  tree  structure  is  used  to  form  the  front  face  list.  The  ordering  of  the 
tree  is  arranged  such  that  the  area  of  the  father  face  is  smaller  than  the  face  of  the 
two  sons.  Figure  13.5  shows  a  binary  tree  of  this  form.  The  area  of  each  face  and  its 
position  are  zilso  shown  in  this  figure.  It  is  noted  that  the  position  of  the  two  sons 
are  located  at  Jsonl  =  I  father  x  2  and  Ison2  =  I  father  x  2  +  1,  respectively.  For 
example.  Face  D  located  in  position  4  has  two  sons  located  at  position  8  (H)  and  9 
(I).  A  face  should  be  added  or  deleted  without  altering  the  bineiry  tree  structure.  The 
ideas  of  the  heap-sort  and  heap-search  algorithms  (Ref.  46)  are  used  in  UNIVG. 

a)  Adding  a  new  face  to  the  face  list: 

A  new  face  is  always  added  first  at  the  end  of  the  tree.  Then,  the  internal  order 
of  the  binary  tree  is  rearranged  by  Somparing  the  face  areas  of  fathers  and  their  sons. 
This  procedure  can  be  described  as  follows: 

1)  A  face  J  with  area  4.0  is  to  be  added  to  a  bineiry  tree  in  Fig.  13.6. 

2)  Place  J  at  the  end  of  heap  list  (10  in  this  case). 

3)  Find  the  father’s  position  of  10  by  using  I  father  =  Isonjl.  Therefore  the 
father’s  position  is  5  and  the  face  is  E  in  this  case. 

4)  Compare  the  area  of  faces  J  and  E. 

5)  If  the  area  of  J  is  less  than  that  of  face  E,  interchange  the  position  of  father 
and  son  (see  Fig.  13.6a). 

6)  Unless  J  is  at  position  1  (top  of  the  list)  go  back  to  (3). 

b)  Delete  a  face  from  the  front  face  list: 

A  face  can  be  removed  from  any  position  from  the  tree.  Then,  the  internEil  order 
of  the  tree  is  re-established  by  comparing  the  area  of  the  face  of  the  father  aind  son. 
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.  -  V-  V,  far<»  is  to  be  removed  from  the  tree.  For 

1)  Find  the  position  from  which  the  face 

we  een  consider  face  A  at  position  1  m  F.g.  13.6b. 


2)  Place  the  face  stored  at  the  end  of  the  tree  at  this  position  (1  for  this  case). 

3)  Find  the  two  sons  location  by  using  fsonl  =  //other  x  2  and  /son2  = 
I  father  x  2  + 1- 

4)  Deteratine  if  the  father  and  son's  portion  should  be  changed  by  checking 
the  area  of  the  father  face  and  two  sons’  faces. 


if 

{Areaifather)  <  miniArea{sonl),Area{son2))  :  no  change 


else 

change  father’s  position  with  the  smaller  area  son’s  position. 


5)  Repeat  step  4  until  no  change  is  reqired  (Fig.  13.7). 

In  this  binary  tree,  the  face  with  smallest  area  will  always  remain  at  the  top  of 
the  list  and  can  be  used  as  the  next  surface  to  be  removed. 

13.1.2  List  for  Points  in  the  Neighborhood 

An  octree  data  structure  is  used  to  efficiently  locate  points  in  a  neighborhood 
The  fist  octant  is  determined  from  the  boundary  of  the  computatmnal  domarn.  At 
r^eight  points  are  stored  in  each  octant.  U  a  ninth  point  falls  into  an  octant,  then 
it  is  s^tbdivided  into  eight  smaller  octants.  This  procedure  is  continued  unt.l  an  oct„ 
.nth  vacant  storage  is  found.  Figure  13.8  iUustrates 

added  and  it  falls  into  octant  1  which  already  contmns  e.ght  pomts  W 

Therefore  octant  1  will  be  divided  into  eight  smaller  octants  (octants  2-9).  Th  y 

“mt;:tholdpointsD.F.Farerelocatedino^^^ 

.V  rOCTRtll  MXOCT)  is  defined  to  store  the  pomts.  The  variable  MX 

n  —  a.  -a  —  .00,  - 
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mformation  is  stored. 


IOCTR{n,IOC)  = -I  :  the  octant  is  full 
lOCTRill,  IOC)  =  0  :  the  octant  is  empty 

IOCTR(ll,  IOC)  >0  :  the  number  of  points  stored  in  this  octant 

IOCTR{l  :  8,  IOC)  :  point  numbers  are  stored  here  if  I0CTR{11,  IOC)  >  0 


The  maximum  and  minimum  x,y,z  are  also  stored  for  each  octant.  By  using  this 
octree,  the  neighboring  points  within  some  specific  distance  of  a  given  point  (i,y,  z) 
ran  easily  be  determined. 

13.1.3  Linked  List  Between  Point  and  Faces 

In  an  imstructured  grid  generation  code,  it  is  important  to  determine  the  faces 
that  surround  a  given  point.  In  order  to  do  that,  an  address  pointer  array 
•^LPION{IPONTy 

for  each  grid  point  IG  is  first  allocated.  The  faces  surrounding  point  “/(?”  are  saved 
in  an  array  LFAPO{8,IG).  The  procedure  to  find  faces  surrounding  a  given  point 
IG  is  described  below; 

1)  Find  the  address  related  to  the  given  point  IG. 

lADRES  =  LPION(IG) 

2)  LFAPO{S,  I  AD  RES)  >  0  defines  number  of  stored  faces. 

LFAPO{S,  lADRES)  <  0  implies  that  next  face  number  surrounding  the  grid 
point  IG  is  continued  at  the  address  abs{LFAPO{S,  lADRES)). 

3)  LFAPO(l  :  7,  lADRES)  =  0  defines  an  empty  location 

LFAPO(l  :  7,  lADRES)  >  0  denotes  face  number  for  face  which  surrounds  the 
point  IG.  This  method  is  illustrated  in  Fig.  13.9. 
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13.2  Check  for  Intersecting  Faces 

In  th.  pr«ess  of  generating  cells,  no  two  faces  mnst  intersect  each  other.  This 

1  -r  of  either  face  intersects  the  other  face  (Fig 

condition  can  be  satisfied  only  .t  no  s.de  o  e^er 

ter#.  13  lO'i  A  total  of  six  conditions  (One  each  t  ,  at 

ILS'  have  to  be  checked  to  make  sure  these  two  faces  do  not  cr^s  each  other, 
tnangles.)  .  .  .  ^  .  detennine  whether  a  side  53  (connecting  X4 

With  the  notation  given  in  the  figur  ,  line  ('iiXs')  with 

a  av  f  n  n  (rt  x->  itV  The  inteTsection  point  oi  line.  (^X4  sJ 

and  X5)  intersects  the  face  gm  (ii 

the  plane  defined  by  ii,  X2  and  X3  is  given  by 

-  .  fqi  X  ga)  -(^1 
=  {h^h)-h 

The  segnrent  (f.fs)  intersects  the  triang^  ‘"f 

itself.  This  procedure  leads  to  three  conditions  for  each  tnangl  . 

are  satisfied,  then  the  two  faces  do  not  cross. 

13.3  Grid  generation  for  viscous  flow 

1)  Follow  the  first  three  steps  of  the  inviscid  grid  generation. 

2)  Generate  the  boundary  of  the  viscous  region  by  moving  the  patch 

its  normal  a  specified  distance  “d”  (Fig.  13.11).  Since  the  normal  at  a  p 
not  uniquely  specified,  an  average  of  the  normals  to  A'  tn-gles  “.at  m^t^_^ 
the  given  point  is  determined  by  taking  the  average  o  e  norm 
adjalt  triangles  that  differ  the  marimum.  A  point  is  then  generated  at  a  given 

distance  along  this  normal. 

3)  Employ  an  elliptic  smoother  to  smooth  the  boundary. 

4)  A  smooth  surface  enclosed  by  this  boundary  is  generated  using  ‘"“f 

'  Lpolation.  This  process  involves  creating  a  map  between  the  local  surface 

coordinates  (s,  t)  and  the  cartesian  coordinates  (x,  y,  z). 
step  4.  This  is  done  by  locating  points  on  the  ne 
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locations  as  the  (s,t)  coordinates  of  the  vertices  of  the  triangles  on  the  body 
surface  patch. 

6)  At  this  stage  the  newly  generated  points  for  all  the  body  surface  patches  are 
checked  for  positive  distance  from  the  body  surface  (Fig.  13.12).  Points  with 
negative  distance  from  the  body  surface  are  corrected  by  a  trial  and  error  pro¬ 
cedure. 

7)  Corresponding  triangles  on  the  two  surfaces  are  connected  to  form  triangular 
prisms.  Any  occurance  of  twisted  prisms  is  rectified  by  moving  the  points  in¬ 
volved  using  an  elliptic  smoother.  The  procedure  for  checking  whether  a  trian¬ 
gular  prism  is  twisted  or  not  is  given  in  Figure  13.13. 

8)  The  triangular  prisms  are  broken  into  smaller  prisms  to  ensure  that  the  flow  in 
the  boundary  layer  region  is  properly  resolved  (Fig.  13.14a). 

9)  The  smaller  tringular  prisms  are  divided  into  tetrahedra  (Fig.  13.14b). 

Figure  13.15  shows  the  body  surface,  smooth  surface  and  projected  triangles  of 
the  space  shuttle  nose  region. 
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Figure  13.1  Surface  grid  generation 
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1)  start  with  (A,B). 

2)  There  are  2  choices: 

a)  A  can  be  connected  to  B’ 

b)  B  be  connected  to  A’ 

3)  compute  lengths  of  AB’  and  BA’  and  choose  the  shorter  one;  for  e,g.,  AB’. 

4)  since  B  is  not  connected  yet,  continue  the  process  with  (A,B). 

5)  2  choices  : 

a)  A  can  be  connected  to  C’ 

b)  B  can  be  connected  to  B’ 

6)  compute  lengths  of  AC’  and  BB'  and  choose  the  shorter  one;  for  c.g.,  BB’. 

7)  since  both  A  and  B  are  now  connected,  continue  the  process  with  (B,C). 


Figure  13.2  Recipe  for  generating  triangular  elements 
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_ .  front  Eurfico  (Evidod  into  tringultr  domonts. 

ABC  fc  CBD  :  two  ot  tho  initial  surface  elements. 

E  &  F  :  newly  generated  nodes. 

ABE  ,  BCE  ,  CAE  :  new  dements  that  replace  ABC. 

CBF  ,  BDF  ,  DCF  :  new  dements  that  replace  CBD. 

BFE  &  FCE  :  new  dements  senerated  by  connectins  two  existing 
“place  the  demenU  BCE  and  CBF. 


Figixre  13.3  Advancing  front  technique 


nodes  E  and  F 
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Figxire  13.5  Binary  tree 
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Figure  13.6  Addition  of  new  face  to  binary  tree 
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FOR  AABC,  (AB  x  AC).  AK  >0  has  to  be  satisfied. 
Same  as  all  other  triangle  ACD,  ADE..etc. 


Figiire  13.12  Criterion  for  accepting  newly  generated  point  A 
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14.0  topology  treatment 

This  section  describes  various  topology  issues  related  to  connectivity  of  nodes 
into  cells,  cells  common  to  a  node,  cells  with  a  common  face,  how  cells  may  be 
subdivided,  how  links  may  be  removed,  etc. 

14.1  NqHp!  Locations  and  Nodes  of  Cell 

The  computational  unstructured  mesh  is  defined  by  1)  specifying  the  set  of  nodes 
and  their  locations  {x,y,z  coordinates  and  optionally  their  nodal  velocities  i,y,i), 
and  2)  specifying  the  set  of  nodes  that  comprise  each  cell.  These  node  numbers 
are  called  the  global  node  numbers.  Each  cell  has  an  ordered  collection  of  nodes 
numbered  locally.  When  the  global  nodes  of  each  cell  are  provided,  they  must  be 
given  in  the  same  order  as  the  local  node  numbers.  Therefore,  given  any  local  node 
number,  the  corresponding  global  node  number  is  easily  determined. 

14.2  Elimination  of  Redundant  Nodes  and  Cells 

Let  us  make  a  set  of  all  nodes  that  are  part  of  at  least  one  cell  in  the  total 
collection  of  cells  given.  This  may  be  a  subset  of  the  collection  of  nodes  for  which 
the  locations  have  been  provided.  The  extra  nodes  in  the  larger  set  are  those  nodes 
not  needed  to  define  any  cell.  These  nodes  (and  their  locations)  may  have  been 
generated  during  the  mesh  generation  process  but  were  not  used  because  their  inclu¬ 
sion  would  have  resulted  in  cells  that  were  somehow  unsatisfactory.  The  ability  to 
remove  redundant  nodes  is  provided  for  in  the  unstructured-grid  UNIVERSE-senes 

formulation. 

During  the  process  of  removing  links,  a  topic  to  be  covered  later  in  this  section, 
some  cells  may  become  degenerate.  This  may  happen  in  two  ways:  1)  the  number  of 
distinct  global  vertex  node  numbers  is  less  than  four;  2)  the  number  of  distinct  vertex 
node  locations  is  less  than  four.  Case  1  deals  with  a  “logically”  degenerate  cell  and 
case  2  with  a  “physically”  degenerate  cell.  Removing  these  cells  become  redundant 
for  the  purpose  of  computing  the  dependent  variables.  Removing  such  redundant 
cells  corresponding  to  case  1  is  available  in  UNIVERSE-senes  unstructured-grid  for¬ 
mulations  as  part  of  various  “grid  editing”  options. 
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1 4.a  Cells  of  Given  Node 

Given  the  nodes  of  each  cell,  the  inverse  map  can  be  constructed.  This  provides 
knowledge  of  which  cells  include  the  given  node.  While  each  cell  has  the  same  number 
of  nodes  (assuming  the  same  cell  type  for  every  cell),  the  number  of  cells  that  include 
a  given  node  varies  with  each  node.  We  have  verified  for  linear  elements  (only  vertex 
nodes)  that  the  sum  of  cells  associated  with  each  node  is  equal  to  the  number  of  cells 
times  the  number  of  nodes  per  cell  which  is  a  somewhat  surprizing  statement  at  first 
glance. 

14.4  Common  Faces 

Each  face  is  made  up  of  a  specific  set  of  local  node  numbers,  each  of  which  is 
mapped  onto  a  global  node  number.  Using  the  knowledge  of  nodes  of  a  cell  and  cells 
of  a  node,  the  cells  that  share  a  common  face  can  eaisily  be  identified. 

Faces  are  numbered  sequentially  with  cell  number.  The  first  six  faces  belong  to 
the  first  cell,  the  next  six  faces  belSiig  to  the  second  cell,  etc.,  for  the  hexahedral  cell 
type.  Therefore,  face  numbers  axe  related  directly  to  cell  numbers.  A  face  shared  by 
two  cells  sports  two  face  numbers,  one  that  identifies  it  as  part  of  the  first  cell  and 
another  that  identifies  it  as  part  of  the  second. 

A  given  face  of  a  given  cell  has  a  number  of  vertex  nodes.  Each  vertex  node  is 
associated  with  many  cells.  Out  of  these,  one  cell  is  the  cell  pointed  to  by  the  face 
number  being  considered.  There  is  not  more  than  one  more  cell  which  is  common  to  all 
nodes  of  the  face.  In  this  fashion,  the  cells  that  share  a  given  face  can  be  determined. 
It  then  becomes  a  relatively  trivial  matter  to  identify  which  face  of  the  second  cell 
is  identical  to  the  face  being  considered.  Thus  the  face  number  corresponding  to  the 
neighbor  cell  can  be  identified,  given  the  face  number  of  a  particular  cell.  Of  course, 
at  a  boundary,  a  face  can  only  be  part  of  one  cell,  the  interior  cell. 

14.5  Cell  Types 


A  major  enhancement  added  during  Phase  III  was  the  ability  to  deal  with  dif¬ 
ferent  cell  types  in  the  mesh.  In  fact  any  cell  could  be  of  any  of  the  types  allowed  by 
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TJNIV  (tetrahedral,  triangular  prism  or  hexahedral).  The  user  specihes  aU  the  types 
that  are  expected  to  be  encountered  in  the  mesh.  At  the  time  of  reading  in  the  mesh. 
UNIV  sets  up  an  cell  type  identifier  for  each  cell.  Each  cell  type  has  an  associated  C 
structure  that  completely  describes  the  local  book-keeping  information  for  that  type. 
Whenever  necessary,  a  cell  refers  to  its  type.  This,  in  turn,  points  to  all  attributes  of 

that  type. 


14.6  Cell  Division 

For  various  reasons,  a  cell  may  have  to  be  divided  into  two  cells.  In  Phase  I, 
cell  division  capability  had  been  developed  for  the  triangular  prism  cell  type  for  use 
in  two-dimensional  inviscid  store-tracking  studies.  Figure  14.1  displays  a  situation 
where  one  face  of  a  triangular  prism  cell  has  been  divided.  If  the  face  is  common  to 
two  cells,  this  leads  to  both  cells  being  subdivided  into  two  cells  each. 

The  following  book-keeping  details  must  be  kept  track  of. 

1.  New  nodes  are  introduced  (with  global  numbers  greater  than  the  existing  maxi¬ 
mum). 

2.  New  cells  must  be  constructed  (with  cell  numbers  assigned  to  be  greater  than 
the  existing  maximum  values). 

3.  Two  existing  cells  must  be  redefined  to  account  for  their  being  made  of  different 
nodes  after  subdivision. 

4.  “Cells  of  node”  must  be  recomputed  in  the  vicinity  of  the  subdivided  cells. 

5.  “Face  to  face”  correspondence  must  also  be  reestablished  locally. 

Such  a  grid  “editing”  capability  has  been  developed.  The  cell  division  capability 
is  currently  not  being  used  (at  end  of  Phase  III). 

14.7  Link  Removal 

It  is  sometimes  desirable  to  remove  links.  In  Fig.  14.2,  the  short  link  ab  is  an 
example.  The  Unk  can  be  removed  by  moving  one  of  its  nodes  towards  the  other.  In 
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this  e^ple,  cells  2  and  5  become  degenerate  and  must  be  removed  from  considera¬ 
tion.  Node  b  merges  with  a  and  therefore  a  is  used  to  replace  b  in  all  cell  definitions. 
Ckmsequently,  node  b  becomes  redimdant  and  can  be  removed  from  consideration  if 
necessary.  The  book-keeping  steps  are  given  below. 

1.  Node  b  for  all  cells  associated  with  b  must  be  replaced  with  node  a  (nodes  of  cells 
3  and  4  must  be  redefined). 

2.  CoUapsed  cells  2  and  5  must  be  removed  from  the  database. 

3.  K  node  b  is  not  removed  from  database  and  must  be  used  for  some  purpose,  the 
location  coordinates  of  b  must  be  modified  to  coincide  with  node  a’s  position. 

4.  “Cells  of  node”  and  “face  of  face”  must  be  recomputed  in  the  local  region. 

Such  a  grid  editing  capability  has  also  been  developed.  The  link  removal  capa¬ 
bility  is  currently  not  being  used  (at  the  end  of  Phase  III). 

1  4.R  Octree  Sort  and  Search 

The  octree  sort  and  search  procedure  applied  to  spatial  data  is  the  three- 
dimensional  equivalent  of  the  binary  sort  and  search  for  one-dimensional  data.  We 
exploit  this  heavily  to  search  efficiently  for  a  point  nearest  to  a  given  point,  a  cell 
that  contains  a  given  point,  intersection  of  an  interior  surface  with  cells,  association 
of  mesh  point  location  with  surface  boundary  location,  etc. 

We  begin  by  sorting  avmlable  mesh  points  into  the  octree  spatial  data  structure. 
Typically,  we  restrict  the  maximum  number  of  physically  distinct  mesh  points  in  any 
given  octree  leaf  to  8.  Duplicate  mesh  points  (different  points  at  the  same  physical 
coordinates)  are  included  and  do  not  have  to  be  discarded  in  our  implementation. 
In  the  following  discussion,  “leaves”  are  terminal  branches  of  the  octree.  Nontrivial 
database  information  is  stored  only  on  octree  leaves.  Octree  branches  that  are  not 
leaves  serve  only  in  the  role  of  a  “container”  pointing  to  the  next  level  of  branching 
or  subdivision.  The  two  basic  procedures  that  can  be  applied  to  the  octree  database 

are  as  follows: 

(a)  Given  a  target  point,  find  the  octree  leaf  that  contains  it. 


166 


NAWCWPNS  TP  8220 


(b)  Given  a  “box”  (of  given  dimensions),  find  all  the  octree  leaves  that 

intersect  with  it. 

1)  By  making  use  of  (a)  and  (b),  it  is  easy  to  determine  the  mesh  point  (that  is  in  the 
octree  database)  nearest  to  a  given  target  point  (that  may  or  may  not  be  in  the  octree 
database).  This  is  used  repeatedly  in  our  grid  generation  method  for  constructing 
unstructured  grids  (with  tetrahedral  or  triangular  prism  cells). 

2)  We  also  construct  an  auxiliary  database  from  the  octree  database.  This  is  the 
list  of  mesh  cells  that  are  associated  with  each  octree  leaf.  In  its  crudest  form,  for 
each  cell  its  “container”  box  is  constructed  and  procedure  (b)  above  is  used  to  find  all 
leaves  that  intersect  with  this  box.  The  cell  under  consideration  is  added  to  the  Hst  of 
“assodated”  cells  for  all  these  leaves.  From  this  information,  the  following  procedure 
ran  be  constructed: 

(c)  Given  a  target  point,  find  the  mesh  cell  that  contains  it. 

3)  This  is  accomplished  by  first  using  (a)  to  determine  the  octree  leaf  that  contains 
the  target  point  and  then  identifying  all  the  cells  that  are  associated  with  this  leaf. 
These  and  only  these  cells  are  sezirched  to  determine  if  the  target  point  is  contained 
within.  If  none  of  the  cells  contains  the  point,  the  “nearest”  cell  from  this  list  can 
be  identified  as  that  cell  for  which  the  normal  distance  from  the  point  to  any  of 
the  cell’s  faces  is  the  shortest.  This  normal  distance  is  defined  to  be  positive  when 
the  point  is  outside  the  cell  and  negative  when  it  is  contained  within  the  cell.  This 
“positive”  distance  will  come  in  handy  in  the  later  section  dealing  with  the  idea  of 
the  “alternate”  boundary  condition. 

4)  Similar  to  the  database  in  Item  2  above,  we  also  construct  (when  necessary)  an 
auxiliary  database  that  lists  the  surface  elements  (or  surface  cells)  that  are  associated 
with  each  octree  leaf.  These  surface  cells  are  used  to  describe  the  geometry  of  physical 
or  fluid-dynamic  entities  that  are  either  at  the  boundary  or  in  the  interior  of  the  flow- 
field  (and  computational)  domain.  They  may  also  be  at  mesh-block  boundaries. 
Using  this  database,  the  following  two  procedures  may  be  defined. 

(d)  Given  a  set  of  surface  elements  (or  even  a  single  one),  find  the  mesh  cells  that  are 
associated  with  the  same  octree  leaves  that  the  surface  element(s)  are  associated 
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with. 

(el  Given  a  set  of  mesh  cells  (or  even  a  single  one),  find  the  surface  ceUs  that  are 
associated  with  the  same  octree  leaves  that  the  mesh  cell(s)  are  associated  wth. 

5)  Procedures  (d)  and  (e)  lead  to  a  very  efficient  way  of  isolating  a  few  surface  elements 
that  correspond  to  a  given  mesh  cell  and  can  be  used  to  define  a  method  to  find  the 
intersection  of  surface  elements  and  mesh  cells: 

(f)  Given  a  mesh  cell,  find  the  intersection  of  the  edges  of  this  cell  and  a  given 
surface. 

6)  Procedure  (f)  can  be  used  to  describe  the  intersection  of  the  mesh  and  the  surface 
geometries.  This  can  be  used  to  generate  interior  or  boundary  condition  specihcations. 
It  can  also  be  used  for  postprocessing  applications  where  it  is  necessary  to  construct 
the  solution  on  an  arbitrary  cutting  plane. 


After  cell  division 


Figure  14.1  Cell  division  for  triangular  prisms 
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Figure  14.2  Link  removal  for  triangular  prisms 
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